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ABSTRACT 


The power spectral methods of generalized harmonic analysis 
are extended to determine the higher order moments above the 
second of the probability frequency functions for either the in- 
put forcing function or the output response of a linear system 
Several techniques are then indicated for determining the corre- 
sponding probability distributions from the calculated moments 
The conditions under which a non-Gaussian forcing function and 
a Gaussian output response can exist are indicated and discussed 
The most significant requirement is that the bandwidth of the 
input power spectrum is large compared with the bandwidth of 


the system transfer function 


INTRODUCTION 


, ] SHE DETERMINATION OF THE RESPONSE of a dynamic 
system due to random disturbances, such as an air- 


plane subjected to atmospheric turbulence, buffeting, 
landing impacts, taxiing over rough surfaces, etc., 
has recently received consideration with the point of 
view of applying the concepts and knowledge of sto- 
chastic processes for describing the statistical char- 


acter of the problem in a more concise mathematical 
form. In references 1, 2, 3, and 4 investigations are 
presented of the application of power-spectral methods 
of generalized harmonic analysis to several of these 


airplane response problems. Although all of the studies 
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realize the necessity of describing the disturbance func 
tion statistically, its description is usually restricted 
to knowledge of only the correlation function or its 
Fourier transform, the power-spectral density func 
If the probability distribution of the disturbance 
the 


tion. 
function is 
power-spectral density approach is sufficient for de 
scribing not only the probability distribution of the 
disturbance function but also the output response func 
For the case of non-Gaussian disturbance func 


Gaussian, then for linear systems 


tion.! 
tions, this information by itself does not, in general, 
adequately describe the statistical properties of either 


the input or the output functions. 


The purpose of the present paper is to present an ex 
tension of the power-spectral density concepts to the 
case of non-Gaussian distributions of the disturbance 
or response functions. In addition, it will be shown 
that, for the case of a lightly damped dynamic system 
and a flat input power spectrum (relative to the transfer 
function of the dynamic system), the input function 
may be non-Gaussian but the response function would 
Limited experimental evalu 
probability dis 


be, in general, Gaussian. 
ation of an airplane's acceleration 
tribution while penetrating atmospheric turbulence’ 
indicated that, the measured cumulative fre 
quency functions were approximately Gaussian, the 
deviations from the for the 
largest accelerations; consequently, proper designa 


while 
Gaussian were largest 


tion of the critical loads experienced by the airplane 
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based on the Gaussian distribution appears question- 
able. In addition, the degree of proximity of the prob- 
ability distribution to Gaussian appeared to be a func- 
tion of the dynamic stability of the airplane, thus indi- 
cating that, in general, the atmospheric turbulence is 
probably non-Gaussian and that the usefulness of the 
Gaussian distribution to approximate the response 
probability distribution would depend upon the transfer 
function of the airplane. Thus, for the general case it 
would appear desirable to assume non-Gaussian fre- 
quency distributions for both the input and the output 
but to realize the usefulness of the Gaussian distribu- 
tion as a good first approximation in some instances. 


SYMBOLS 


t = time, sec 
ty, ts = dummy variables of integration 
T1, T2 = arbitrary displacements in time, sec 
ii = arbitrary value of /, sec 
x(t) = random input forcing function 
y(t) = random output response function 
h(t) = response of linear system to a unit impulse 
forcing function 
f = frequency, cycles per sec 
H(f) = transfer function or frequency response of 
%~ 
linear system, H(f) = | h(t)e —I-*"" dt 
J0 
Gf, fe, f,) = power spectrum of input generalized to 
multiple frequencies 
GA fis Ja f,) = power spectrum of output generalized to 
multiple frequencies 
g(x) = Gaussian probability frequency function 
(x) = Gaussian probability distribution func 
x 
tion, P(x) = f g(x) dx 
f(x) = non-Gaussian probability frequency func 
tion 
F(x) = non-Gaussian probability — distribution 
ad 
x 
function, F(x) = | f(x) dx 
My = mth central moment of probability fre 
quency function 
m = mean of probability frequency function 
o = standard deviation of probability fre 
quency function 
E( ) = mean value of a function of a random vari 
able 
=V-1 
y(t) = average value of y(t), y(t) = 
. l 7 
lim — y(t) dt 
T+ 2TJ-1 
"3 = third absolute moment about the mean 
of probability frequency function 
R( 71, 72, T,) = correlation function generalized to mul 
tiple arbitrary displacements 
a = arbitrary constant 
C; = defined by Eq. (A-3) 


Ay = defined by Eq, (A-7) or (A-8) 
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ANALYSIS 


Relations Between the Input-Output Moments of the 
Probability Distributions for a Known Transfer Function 
Consider the problem of determining the moments 

of the probability distributions for either the random 

input, x(t), or output, y(/), given one or the other and 
the transfer function of the system, //(f). If the ran- 
dom input is stationary over sufficiently large time in- 
tervals and the system is linear, the ergodic hypothesis 
becomes valid, implying that time averages are equiva- 
lent to ensemble averages for both the input and out- 
put. This feature of the random process can simplify 
the calculation of the moments of the probability dis- 
tribution through the use of the correlation function. 
This function, together with its Fourier transform, the 
power-spectral density function, provides a technique 
of calculating the second moment of either the input 
or output given the transfer function of the system 
and the power-spectral density of either the output or 
input functions. Denoting, respectively, G,(f) and 

G,(f) as the power-spectral density functions of the 

input and output, a useful relation for determining the 


desired second moments is 
G,(f) = G,(f) HT(f) - (| 

where the bars indicate the modulus of the function. 
A derivation of Eq. (1) is given in reference 5. The 
input or output spectrum can be obtained algebraically 
from Eq. (1) relatively simply; the corresponding second 
central moments or standard deviations, o, and o,, are 
determined readily from the integral of the power 


S ieee 
be, =— e." = ;f G,(f) df 


l oe : 
jin, = oy = / Gety) us| 
uy y 9 ‘ . 
9) 


but since G,(f) and G,(f) are even functions, Eqs. (2 


spectrum. 


simplify to 


. 


Mey, = Sy" = / G,(f) us| 


If the input and output probability functions are Gaus 
sian——that is, in a form given by Eu. (4), 


g(x) = (1/0, V 2x) e~ 242" (4) 
then all that is required for determining the probability 
frequency function, g(x), is the standard deviation, 
O;. 

Examination of Eq. (1) indicates that the effects of 
phasing in the transfer function are obviously not re 
quired in determining the probability distributions for 
the Gaussian that is, 


case. For the general case 
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the non-Gaussian case—it is known that the higher 
order moments (third, fourth, fifth, etc.) become sig 
nificant in determining the probability distributions, 
thus implying a need for the effects of the transfer 
function phasing. It would be desirable then to de 
termine the equivalent relation given by Eq. (1) which 
will allow for the calculation of the higher order mo 
ments, as well as properly accounting for the effects of 
these 


transfer function phasing. A development of 


relationships is shown as follows: 


Consider the following convolition integral relating 
the input and output responses in the time plane 


7 


y(t) = | aT — t))h(t,) dt, » 
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where /i(t) is the response to a unit impulse and /(f 

Qfort < 0. 
placement, 7, of the input response is 


The output response for an arbitrary dis 


vit + 7m) = / v(t + 7; to)h(ts) dt 6 


- 


Similarly, for a second arbitrary displacement, 7, 10 


the input response, the output response becoines 


v(t + 7) = / V(t + Te ta)h(ts) dt 7 


Since the ergodic hypothesis is assumed valid, a three 
dimensional correlation function can be formed for the 
input function and output response from Eqs. (5), 


(6), and (7). For the output, this correlation function 


is 
R,(11, t2) = v(t) v(t + ryt + 72) (S 
where the bar denotes an averaging process, as follows: 
ee 
vit)v(t + r)¥(t + 72) = lim 1 / v(t)yv(t + ri)v(t + re) dt 
roo 2] / 


or substituting Eqs. (5), (6), and (7) into Eqs. (S) gives 


R,(n, T2) = | / / x(t — txt + ory — be)x(t + re — ty) ACA te) A(ts) dt, dts dt v 


where the input correlation function is 


RwA(7, — te +h, t2 — t 


+ fy) =a gif 


hox(t + lo)x(t + To f (10 


Substitution of Eq. (10) into Eq. (9) gives the important relationship for R,(71, T2), 


R, (11, T2 = | | | R,(7 _— ts + fi, >. ts + tL )A(t))h(te)ht \dtydtedt 11 


As is shown in the literature for the two-dimensional case,> a Fourier transform pair also exists for the power 


spectral density and correlation functions for the three-dimensional case. 


These relations are 


G(fi, fe) = 2 f / R (TT), T2)e PaeATetE Tae drdr» (12 


1 [{ oar ~~ 
R(n, yy = / / G( fi, foje~ 
9 


A derivation of Eqs. (12) and (13) is given in Appendix (A). 


of the output response results in the following expression : 


G,| fi, fo) = 2 | / / / / R,(7 =~ be a hh, % — ts 


Substitution of the following variables in Eq. (14) 
Pi = 71 le thy, 


gives the result for G,( fi, / 


df, df. (13 


Substitution of Eq. (11) into Eq. (12) for the case 


+ fiJe h(t )A(ts)A(ts)dtydtodtsd rd re 


(14 
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d pid po: / ” io TI) h(t) dty- | 


or 


Gif, fe) = 2 f / ki he 


but 


Gl fi, fo) = 2 / / R,(pi, Pp») e j2x( fiprt ») dp d po 


(1S) 

H*(fi + fe) = [ Cr hy) iat (19) 
H(fi) = / e2*!2 (to) dte (20) 
(21) 


Hf.) = / e 2? "#5 hits) dts 


where the asterisk denotes the complex conjugate. 
Substitution of Eqs. (18), (19), (20), and (21) into 

Eq. (17) gives the important relationship, 

(22) 


Gif fe) = Gf fo *(fi + fof) HCfs) 


from which the third central moments for either the 
input or output can be determined [see Eq. (13) when 
n= = 0 
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N(ty)A(te)h(ts)dtydtedtsd pid po 16 


| eo" hits) dt f gr 


* h(ts) ds) 


Gifu fe) df dfr 


tg" 
us, = R,(0, 0) = / / 
May = R,(0, 0) = G, | fi, fo) dfidfs (23b 


Note that Eq. (22), which is analogous to Eq. (1) for 
determining the second moment, introduces the effects 
of phasing of the transfer function. 

The analysis shown for the third moment can be ex- 
tended to the case of the fourth moment. Since the 
technique involved is essentially the same, only the 
significant results are given. The correlation func- 
tions for the input and output are defined by the fol- 
lowing relations: 

R,(71, 72, T2) = X(t) X(t + 71) x(t # re)x(t + 73) / (09 
R,(71, t2, 73) = v(bw(t + r)y(t + ro)v(t + 73) \ - 
The Fourier transform pair relating the correlation func 
tion and the power spectrum for either the input or 
output is 


G(fi, fe, fs) = 2 | / / R(7,, Te, T3)e $2 a fala tfaret fora ‘dtd tod 73 (25 
| ee i " _ oo | | 
R(n, T2, T3) = * / / / Gf 1, fa, fade j2 x (fir, + fore "df d fod fs (26 


The input-output power-spectrum relation, which corresponds to Eqs. (1) and (22), for determining the fourth 


moment is 


Gy (fir fr, fs) = G,(hfi, fo, As) ET* (fi + fo + fs )ET( fi )IT( fo) (fs) (27 


The fourth central moments for the input and output are 


pa, = R,(0, 0, 0) = - / / / Gi( iy fo, fs) df dfodf (28 
us, = R,(0, 0, 0) = if / | Gy( fi, fo, fs) dfidfodfs (29 


From an inspection of the results for the third and fourth central moments, extension of the results to higher 


order moments becomes obvious; the respective generalized equations for Eqs. (24), (25), (26), (27), (2S), and 


(29) are 


R,(n, ee 


Ryn, Ti + os 


T,) = x(t)x(t + r)x(t + 7)... 


Tr) = y(t)y(t + ri)y(t + 72)... 


~ Xt + t,) (30) 


y(t + 1,) (31) 
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G(ft f oP = 2f | =F Rr, , ee T.3 . 7 dr, dr dr 4 
R Ti,» T2, Ty = se 6 G hei ° oe is adf\da} day} ts) 
G, (fi, fe Fr) = CAF. 7 f, II*( fi + fo 4+ f, ET f,)LT( fe TT(f 4 
: ° ‘ ; 
M = R,(0, 0 Q / | / G f f f dt, df if 5) 
) 
; | ; ; 0 +) dtd 
) = f f d}hidfo, . a in 
mM °¢ (), 0 : I WidJe 
indicates that An excellent dissertation of these techniques, as well 


Inspection of Eqs. (30) through (36 


necessary for determining moments 


fourth become extremely laborious. 


the calculations 
higher than the 
Consequently, a method for determining the probability 
distribution that requires a minimum number of mo- 
ments is desirable. 

OF NON-GAUSSIAN PROBABILITY 
MOMENTS 


REPRESENTATION 
DISTRIBUTIONS IN TERMS OI 


One of the well-known techniques of estimating a 
probability distribution is in terms of its moments. 
Chis method, introduced by Pearson® * consists in de 
termining a number of moments that are functions of 
unknown parameters describing the probability func 
tions. For an equal number of known moments and 
parameters to be estimated, a set of simultaneous equa 


tions mav be formulated. Solution of the equations 


with respect to the unknown parameters provides a 
means of expressing an approximation of the probabil 
itv frequency function in terms of the measured mo 
ments. More refined methods for developing expres 


sions for the probability functions in terms of moments, 


as well as the implications involved, are available. 

; 1 ps d*o(x 1 fms dig(x) 10 

Io) = 2) - 5 .”CtC — 3 se 
ao” dx t \ G4 dx* 8) 


where the terms on any one line are of the same order of 


magnitude. The corresponding expansion for F(x) in 
terms of &(x) can be obtained from Eq. (37) by replac 
ing g(x) by P(x For a more detailed discussion of 


the Edgeworth series, the reader is referred to the ex- 
cellent dissertation given by Cramér.’”° Examination 


f indicates that, if the second, third, and 


ot Eq. (37 


as a detailed list of references, is covered by Shohat 
and Tamarkin 

One special technique of describing the probability 
distribution and frequency functions in terms of the 
moments is the rhis_ technique 


may be the most practical method for the gust prob 


Edgeworth series 
lem, since probability functions of measured load re 
sponses appear to be approximately Gaussian. A 
series representation of the actual distribution in terms 
of a Gaussian having the same second moment would be 
minimum of terms in the 
Essentially, the Edgeworth 


advantageous, since only a 
series would be required. 

series consists of an expansion in terms of the Gaussian 
distribution, its derivatives, and certain constants. 
These constants are evaluated with the aid of an orthog 
Hermite polynomials in terms of 


onal expansion of 
the standard devi 


moments, u,, and 


and (x), the Gaussian fre 


known central 
ation, o. If we define g(x 
quency and distribution functions, respectively, and 
f(x) and F(x), the observed frequency and distributions 
functions, respectively, then the Edgeworth series can 


be written as 


vi > d' g(x 
S) dx* i 

35mm | Jy(x) 280 (pu 
7 :( _ * dx’ v & 


fourth moments are known, a three-term approximation 
is obtained and appears sufficient for representing 
either the frequency or distribution functions for ap 
proximately Gaussian functions. 
Determination of either f(x) or F(x) by Eq. (37 
together with the technique developed in the previous 


section for determining either the input or output mo 
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ments, allows for an accurate description of the statis- 


tical character of many dynamic problems in a concise 


mathematical form. 
problem of establishing, as well as applying, a rational 
design criterion for the load factor experienced by an 
Old 


designs that are critical for gust loads and have logged 


airplane encountering atmospheric turbulence. 


long satisfactory gust experiences can be used to de 
statistical characteristics of turbulence, 


or ‘‘level’’ for new 


termine the 
as well as to establish a ‘‘reference”’ 
designs that is, to establish the design probability of 
exceeding the design gust load related to past safe ex- 
perience. Analyses of a large number of response time 
histories of accelerations or stresses measured on these 
airplanes while penetrating turbulence can be made in 
the manner suggested herein. The moments and their 
associated power spectrums, as well as the probability 
functions, f(x) and F(x), can be obtained for a variety 
of atmospheric conditions— that is, turbulence over flat 
land, sea, mountains, etc., as well as turbulence at 


various altitude levels. 


The level of probability of exceeding the design gust 
load factor must be selected from previous designs 
that were critical for gust loads. This procedure con- 
sists of calculating the probability function, F(y), from 
records of a large number of airplane responses due to 
flight in turbulence from the techniques presented 
herein. For each design, the desired level of probabil- 
ity that is, the value of F(y) corresponding to Vgesien 
or F(Vgesien) can be determined. If desired, an average 
value Of F(Vgesign) Obtained from several airplanes 
may be chosen as a rational level of probability. Once 
F(Yaesign) 18 known, as well as the statistical properties 
of turbulence, the value of Vgesign for a new airplane 
can be determined from its probability function, F(y). 
The technique of determining F(y) is essentially iden- 
tical to the inverse problem of determining the F(x) 
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function for atmospheric turbulence from measured 
airplane responses. 

INPUT AND 


THE SPECIAL CASE OF A NON-GAUSSIAN 


A GAUSSIAN OUTPUT 


If the bandwidth of the random input, x(f), is much 
wider than the bandwidth of the system transfer func. 
tion, then the output, y(‘), is offen Gaussian. Proof of 
this result is possible with the aid of a theorem due to 
Liapounoff given in Cramér (reference 10, page 215 
Essentially, this theorem states that, if the output is 
a result of a sum of independent variables, 


y(t) = yilt) + volt)... Vat) 


with mean m,, and standard deviation, o,, and the 


vy, 


third absolute moment of y, (¢), about its mean, 


n,,) = E(\y, — my,}*) 
is finite for every vy = 1, 2,3. ...m, where 
eS =4, + &, + wm,” +. > He 
then the sum y(t) = yi(t) + yo(t) +... ¥,(f) is asymp 


totically Gaussian if 


lim 7, o, > O 38 
Assume the following conditions: 
(1) Mean of the input, 7%, equal to zero (no loss of 
generality results from choosing zero mean 
2) The input power spectrum (two-dimensional), 


G,(f), is constant with respect to the system transfer 
function, //(f). [See Sketch (a) ]. 

Let fy be the natural frequency of //(f) and let f, ap 
proach infinity or in the time plane Ar = |. 2f, approach 
zero. For this specified input power spectrum, G,(/ 
it can be shown’ that the standard deviation of the in- 


put is 


> é ” 
Cr” = Wo 2Ar (oD) 


The standard deviation of the output can be expressed 
in terms of o, and the system response to an impulse, 
h(t), with the aid of the convolution integral. Writing 
this integral in finite form, the output standard devi 
ation becomes 


oO, = Yo, h?(t — 7) (Ar)? 10) 
0 
Substitution of Eq. (39) into Eq. (40) gives 


o, = 0S h(n) Ar t] 


2 0 


For the third absolute moment the convolution integral 
can again be used to give the relation 


nm? = >> nth3(t — 7) (Ar)? (42 


0 


If the condition given by Eq. (38) is to be satisfied as 
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approaches infinity (or Af approaches zero), then n, 
must approach zero, since ¢, must always be greater 
than zero |/7(f) is not an impulse]. In order that 
7, approach zero, it is seen from Eq. (42) that ,* must 
ypproach | (Aft)* as Ar approaches zero for finite h(t). 
For the output, y(t), to be Gaussian, n,/0, approaches 
zero if n;° approaches infinity slower than | (Ar)* ap- 
proaches infinity. But is proportional to 


Ar [Eq. (2)], then n, 


since a," 
must approach infinity slower 
than ¢o,* approaches infinity. Thus if the input is 
non-Gaussian but its spectrum is wide compared with 
the spectrum of the system transfer function, then the 
output is Gaussian if the input second moment and 


third absolute moment are such that 


) > O(n, (43 


O(c, 


where the order 0 means the rate as o, or n, converges 
to infinity. 

Consideration of the significance of Eq. (43) indi 
cates that, unless the input probability function con- 
tains an abnormally large amount of skewness, the 
output distribution is Gaussian. Although in many 
instances little or no information may be available on 
the second and third absolute input moments, calcula 
tions for these quantities for satisfying the requirement 
given by Eq. (43) can be made from measured output 
responses employing the techniques developed in the 
present paper. 

For turbulence problems, it is doubtful that its prob- 
ibility frequency function contains sufficient skewness 
to violate the requirement given by Eq. (43). In order 
to obtain some physical insight as to the meaning of 
Eq. (43), an example is shown which will indicate the 
relative order of skewness and standard deviation for 
a well-known unsymmetrical probability function, the 


Rayleigh, defined as follows: 


(for x < 0) | 
(for x > O) 


f(x) = 0 (44 


f(x) = (x/a*)e 


A plot of Eq. (44) for a typical value of a is given in 


Sketch b. 
The standard deviation of f(x) defined by Eq. (44) is 
or =~ av? — (r/2) (45) 
while its mean is 
mM, =av x/2 (46) 


Substitution of the change in variable 


in Eq. (44) gives the following result : 


ae gs eT a 
iu) = = _ ~ ee OE UB 


Gg 


— m,) (47) 


The third absolute moment cubed of f(x) is 
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f(x) 
\ 
oO x 
sKetcH (6b) 
; uM, , 
= ju|° ' e - tS 
) a 


Let v = ua and substitute Eq. (46) into Eq. (48 


the resulting equation becomes 


; T 
S.” = q” 0 9 abe ? \ - 2d (49 
< “ \ 9 é “s 
V x/2 ~. 


Since the value of the integral in Eq. (49) is not a fune 


tion of a, then ,* is proportional to a* or 7,° is propor 
tional to a*. This result, together with the propor 


tionality of o, to a given by Eq. (45), gives the rela 


tionship of the order of o, to n, as 


50 


O(a,) = O(n, ( 


Thus the O(c, *) must be greater than the O(n,), and 
Eq. (43) is satisfied; consequently, a Rayleigh prob 
ability frequency function input that has a wide band 
width power spectrum compared with the system trans 


fer function leads to a Gaussian output 


In the application of the results presented herein for 
turbulence problems, the bandwidth of the input power 
spectrum will be large compared with the system trans 
fer function if the system is only slightly damped. 
Such may be the case of an airplane component such 
as a tail buffeting near its flutter speed or if the air 
plane, while traversing atmosphere turbulence, experi 
ences a large rearward center of gravity movement 
that substantially reduces the damping of its short 
period oscillation. For these cases, experimental in 
vestigation of the character of the second moment and 
third absolute moment of the input probability func 
tions appears desirable in view of the possible simpli 
fications in describing the output probability functions 


if the input is non-Gaussian. 
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In the investigation of deducing the input charac 
teristics from a measured output, it appears advisable 
that, in order to obtain the best estimate of the input 
characteristics, the function should 
have a bandwidth comparable with that of the input. 


If the bandwidth of the system transfer function is ex 


system transfer 


tremely narrow compared with the input, then the ac 
curacy of the input power spectrums would be ques 
tionable at all frequencies except the natural frequency 
additional 


of the function; as an 


system transfer 
all the 


function 


consequence, calculated moments of the in 


would also be question 


put probability 


able. 
CONCLUDING REMARKS 


The higher order moments above the second of the 
probability functions for either the input forcing fune 
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tion or the output response of a linear system are deter 
mined from an extension of power-spectral methods ; 
generalized harmonic analysis. Examination of th 
generalized expressions for determining any order mo 
ment indicates that the calculations necessary for de 
termining moments higher than the fourth become ex 
tremely laborious. Methods are discussed tor describ 
ing probability functions in terms of the moments with 
special interest to one that requires a minimum nun 
ber. For the special case of a non-Gaussian forcing 
function that has a wide bandwidth compared with 
the bandwidth of the system transfer function, it is 
shown that for most cases the output response will be 
The whether the 
output response is Gaussian depends on the relativ 
the third absolut 


Gaussian. final determination of 


order of the second moment and 


moment of the input probability functions. 


APPENDIX (A 


Proof of Fourier transform pair: 


G(fi, fe) = 2 / / R(71, To) € Set) Poul s (A-la 
l ’ P j2 . ‘ 
R(7m, T2) = ; G(fi, foyer” dfidf (A-1b 
Let 
R(7, 72) = lim Cr(11, 72) A-2 
where 
es | 
Crit, T2) oT yr(t)yr(t + ri)yr(t + re) dt A-3 
and where 
Vy v(t (for -7 <t< T) 
vr = 0 (elsewhere) 
The Fourier transform of Cy can be written as follows: 
. : . 2 | . . | 
Cr(m1, T2)e 7” dr,\dtT. = oT c  Vr(t)dt- 
j oes yr(t + 71) dr: / e~ vr(t + raids (A-4 


Let s) = ¢-+ 7, and so = ¢ + 7, and Eq. (A-4) then 


. » 


/ / Cr(n, te)e 2? dridte = ae 


but 


yr(tje~" 


becomes 


‘dt- 


> 


/ gon sears ¥r(S1) dsy: / ee  Yr\S2 dso (A-5 





the: 


Sub 


Fou 


Yor 





re cle ter 
thods ¢ 
Ol the 
der mo 
for de 
DINe ex 
describ 
its with 
n nun 
forcing 
d with 
n, it is 
will be 
ler the 
elative 


bs ute 


of the Aeronautical Sciences, Vol. 20, No. 5, pp. 317-329, May, 


POWER SPECTRAL METHODS OF HARMONIC ANALYSIS 153 

17 (f / Vr(s,)e “ ds, A-7 

! / yr( see 72" ds \-S 

where the asterisk denotes the complex conjugate. Substitution of Eqs. (A-6), (A-7), and \-S) into | q \-5 


results in the following expression for the Fourier transform of Cr (7), 72 


. ' ; 
| / Cr(1m, t2)e 7° dridr 27 Ar*( fi, + fo)Arl fy Arlt A-9 


If we define the power spectrum function as 


G( fi, fo lin (1 T)Ar*(fy + fAr(fiArths 


then Eq. (A-9) becomes 


G(fi, fe ‘ . x ( firitfer 
= = lim Cr(71, T2)e 7 dr ,d ts A-10 


Substitution of Eq. (A-2) into Eq. (A-10) gives Eq. (A-la). Eq. (A-lb) is a direct consequence of the inverse 


Fourier transform and may easily be verified. 
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Effect of Diffusion Fields on the Laminar 
Boundary Laver 


JOHN W. SMITH* 
Bell Aircraft Corporation 


ABSTRACT 


A theory is developed which describes the effect of a general 
diffusion field on the dynamic and thermal characteristics of a 
laminar boundary layer on a flat plate in steady compressible 
flow. Fluid properties are considered as functions of temper 
ature and local concentration of the foreign gas. The diffusion 
field is described by a differential equation that relates convec 
tive and diffusion transfer and which considers diffusion currents 
arising from gradients of concentration and temperature. By 
means of the usual transformations the system is reduced to a set 
of ordinary differential equations, which in turn are transformed 
into a set of integral equations. The latter is amenable to solu 
tion by the method of successive approximations 

The theory and results have bearing on the problem of control 
and reduction of aerodynamic heating at hypersonic speeds 
The special feature of this approach lies in the utilization of dif 
fusion fields for the purpose of reducing the detrimental effects 
of viscous dissipation. Although the theory is adapted to a 
fuller investigation of this problem, the numerical examples 
considered involve mainly diffusion fields of helium, with which 
good results have been achieved at Mach Numbers 8 and 12 
Whereas at the higher Mach Number the influx of heat was prac 
tically eliminated, a reversal in the direction of heat flow has 
been effected at the lower Mach Number 


List OF SYMBOLS 


Cp = specific heat at constant pressurc 

Cj; = constants of integration in Eq. (2.8 

pies 

C; = average peculiar velocity (cf. reference 5 

D\». = coefficient of mutual diffusion 

f(y) = dimensionless stream function [ef. Eq. (3.10) 

F, = functions defined by Eq. (3.18) 

phe 

J, = diffusion current 

k = thermal conductivity 

K, = constants of integration 

J = Mach Number 

n = total particle density 

n, = foreign particle density 

nin = (nj/n) 

p = pressure 

P, = Prandtl Number 

t = (T/T») 

7 = absolute temperature 

“u,v = boundary-layer velocity components in x, y directions, 
respectively 

Xx, ¥ = coordinates parallel and perpendicular to flat plate 
surface, respectively 

a = thermal diffusion constant 


8, = constants defined by Eq. (2.19) 


= ratio of specific heats 
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n = independent variable introduced by Eqs. (2.10) and 
(3.10) 

Xr = { dD, /v 

m = coefficient of viscosity 

v = kinematic viscosity 

p = density 

gn = functions defined by Eq. (3.11 

gni = functions defined by Eq. (3.23 

y = stream function 

Vv, = functions defined by Eq. (3.18 


Subscript 0 generally refers to free-stream quantities 
Subscript aw refers to adiabatic wall conditions 


Primes indicate differentiation with respect to 7 


(1) INTRODUCTION 


£ os BASIC MATERIAL UNIT considered in standard 
theories of fluid mechanics is the fluid particle, 
which is conceived as a volume of fluid having the 
following properties: It is arbitrarily small compared 
with the macroscopic dimensions of the physical system, 
vet it encloses a sufficiently great number of molecules 
so that 
density can be neglected. 
sidered to move with a velocity that at each point is 
velocity 


fluctuations in pressure, temperature, and 


The fluid particle is con- 
equal to the macroscopic flow described 
by the standard equations of fluid mechanics, and its 
trajectories are just the stream lines of laminar flow. 
Now let us suppose that the molecules constituting a 
fluid particle at a given instant are labeled in some way. 
One would then observe a diffusion current of these 
molecules flowing out of the fluid particle in all direc- 
tions, until the given set of molecules has distributed 
itself throughout the entire system. This diffusion 
current is itself a macroscopic phenomenon, observable 
through ordinary instruments, and must not be con- 
fused with the Brownian motion executed by individual 
molecules. The reason that diffusion currents can be 
ignored in most discussions of fluid mechanics arises 
from the circumstance that the constituent molecules 
of the physical system are assumed to be dynamically 
equivalent, in which event the diffusion transfer does 
not affect the characteristics of the flow. 

The local concentration of special molecules, de- 
fined either as a particle density or as a mole fraction, 
may be considered as a parameter defined at each point 
in the flow and constituting a field. This diffusion 
field can be determined theoretically by integration of a 
differential equation that is analogous to the continuity 
equation of ordinary fluid mechanics and which ex- 
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presses the conservation of special molecules. As was 
pointed out above, the special case of diffusion transfer 
involving dynamically equivalent molecules leads to 
great simplifications because the velocity, temper 
ature, and pressure fields are unaffected by the diffusion 
transfer and can be calculated without recourse to the 
diffusion equation. Diffusion of active nitrogen into 
inert nitrogen and mutual diffusion of isotopes of a 
given element are instances of diffusion problems that 
may be considered in this fashion, since the condition 
of dynamic equivalence is closely satisfied. In general, 
however, such parameters as density, specific heat, vis- 
cosity, ete., which enter into the ordinary flow equa- 
tions, depend strongly on the concentration of special 
molecules, for which reason the ordinary flow equations 
cannot be solved without recourse to the equation of 


diffusion. 


The present paper concerns itself with steady-state 
diffusion fields of a foreign gas in the laminar boundary 
laver on a flat plate and its effect on the dynamic and 
thermal boundary-layer profiles. Apart from the 
theoretical interest attaching to the given problem, the 
theory has important bearing on boundary-layer con- 
trol through fluid injection, a topic that has received 
considerable attention in literature.' Most 
theoretical treatments of fluid injection assume dy 
namically equivalent fluids and study effects that are 
entirely different from the effect considered in this 
The now classical report of Schlichting and 
main 


recent 


paper. 
Bussmann,” for instance, studies the effect of 
taining a nonvanishing normal velocity at the wall on 
the velocity profiles of an incompressible boundary 
layer. Other reports have studied the cooling of a 
compressible boundary layer through injection of liquid 
air, which, in addition to the effect considered by 
Schlichting and Bussmann, involves energy transfer 
between the boundary layer and coolant. Although 
the present theory is fundamentally adequate to ac 
count for these various conditions, the sample problems 
given here have isolated for consideration the special 
effect of the diffusion field of a foreign gas on the ve 
locity and temperature profiles of a laminar boundary 
laver. 


This effect appears when the physical characteristics 
of the foreign gas, such as density, viscosity, heat con 
ductivity, and specific heat, differ appreciably from 
those of the primary gas. The qualitative and quan- 
titative character of the effect depends, of course, on 
the particular foreign gas selected. A thorough inves 
tigation of this problem would utilize hypothetical 
gases for which the combination of physical character- 
istics can be varied at will, and the influence of these 
characteristics on the boundary layer profiles could be 
studied systematically. In this fashion one could find 
criteria for the selection of particular foreign gases 
boundary- 


best suited for any given purpose of 


layer control. Because of the rather time-consuming 
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nature of the numerical calculations, especially for the 
compressible case, we have adopted the simpler pro 
gram of selecting for consideration one or two special 
gases that have interesting and possibly useful physical 
characteristics and for which empirical determinations 
of the relevant physical constants have been made. 
Hydrogen and helium, both notable for their low densi 
ties, high heat conductivities, and high specific heats, 
were considered in the isothermal (low-speed) problem, 
and helium alone was considered in the case of com 
pressible (hypersonic) flow. As will be seen in Sections 
3) and (4), helium gas holds great promise as a coolant 
to be applied in flight régimes where aerodynamic heat 


ing would become important. 


The simpler case of isothermal boundary-layer flow 
has been investigated first, and the theory was later 
extended to the more complicated case of compressible 
flow with dissipation of kinetic energy. It was thought 
expedient, for reasons of clarity, to follow the same de 
velopment in the present paper. Section (2), there 
fore, contains a brief account of the simpler problem, 


while the main theory is presented in Section (3 
(2) ISOTHERMAL BOUNDARY LAYERS 


(a) Basic Equations 


The laminar dynamic boundary layer on a flat plate 
is described, as usual, by three differential equations, 
expressing conservation of mass, conservation of A 


momentum, and conservation of y-momentum, respec 


tively, 
O Ox) (pu) + (O Ov) (pz 0 2 
Ou Ou Oo —) 
pu + pi - (« 5 ai 
OX Ov Ov OV 
Op Ov 0 > 3 


Because of Eq. (2.3), the pressure must be constant 


throughout the field. The condition that the boundary 
layer be isothermal implies, therefore, that the density 
p and the coefficient of viscosity uw are functions of the 
concentration parameter only. Suppose 7; denotes the 
particle density of foreign gas, and m denotes the total 
particle density. According to Avogadro's law, m is 
constant throughout the boundary layer, and we can 
dimensionalize m, by referring it to nm. The quantity 
Nw = n, n, which is the local mole fraction of foreign 
gas, will be considered as the basic concentration param 


eter that defines the diffusion field. 


The diffusion current /; is defined as the total quan 
tity of special molecules (expressed here in terms of 
mole fractions) flowing through unit area perpendicular 
to the flow direction in unit time. At constant tem 
perature and pressure and in the absence of body forces, 
the diffusion current can be related to the diffusion field 


through the equation, 
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= 


J, = —Dy grad nyo (2.4) 


where the coefficient of mutual diffusion Dj,» is largely 
independent of concentration. The differential equa- 
tion expressing the conservation of special molecules is 
derived by balancing the convection and diffusion cur- 
rents, 


u(Ono/OX) + v(ON/Oy) = — div J (2.5) 


Using Eq. (2.4) and assuming that 0°/0x? < 0° Oy’, 
the diffusion equation for the isothermal boundary 
layer becomes 

O70 


= Dy : (2 6 
Ov" 


On 0 On 
u v 
Ox Oy 


(b) Self-Diffusion 

The somewhat trivial case of self-diffusion deserves 
brief mention, because it is easily dealt with and has an 
interesting application. If the special molecules are 
dynamically equivalent to the primary molecules, then 
p and uw are constant, and the x-momentum equation 
can be written in the form 


u(Ou/Ox) + v(Ou/OV) = v(O07u/Oy") (2.7) 


where vy is the kinematic viscosity. The coeflicient of 
self-diffusion D,,; has the dimensions of kinematic vis 
cosity and, in fact, has roughly the same magnitude. 
To a first approximation, at least, D,, can be identified 
It is evident, then, from a comparison of Eqs. 

2.7), that the solution for the concentration 


with »v. 
(2.6) and 
profile is obtained as a linear function of the velocity 
profile, 

No = Ci + Cu (2-8) 


C, and C, being constants of integration. 

The fact that the concentration profile for  self- 
diffusion coincides, roughly, with the velocity profile 
for isotherinal flow over a flat plate could be utilized for 
the purpose The 
obvious experiment that is suggested by this result con- 


of boundary-layer visualization. 
sists in the use of active nitrogen to be introduced into 
an inert boundary layer through some permeable wall. 
It is supposed that this can be done in such a manner 
that the boundary-layer flow is virtually undisturbed. 
One can expect that an orange glow will appear in the 
boundary-laver region, the thickness of which is com- 
mensurate with the thickness of the dynamic boundary 
layer, and that the transition point should be con- 


spicuous. 


c) Mutual Diffusion: Method of Solution 


For mutual diffusion, the problem is complicated by 
the fact that the local density and viscosity become 
functions of the concentration parameter my. An im- 
portant simplification is achieved by assuming that 
both # and 7, are functions of the similarity parameter 
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n alone, where y ~ yx This assumption is just; 
fied so long as the boundary conditions permit all rele 
vant quantities to remain constant at y = 0, with th 
exception of the normal velocity v. The latter may 
in general, be permitted to vary inversely as the square 
root of x. When these conditions are satisfied, it js 
possible to reduce the given system of partial differ 
ential equations to a system of total differential equa- 
tions, which is more amenable to numerical integration, 

Let a stream function y, a similarity variable n, and 
a function f(y) be introduced through the following rela- 


tions: 
py OY po OY 
“= : v= — 2.9 
p Ov p ON 
n = (1/2)VV potto, wor 2.10 
Y= V wxmy f(n (2.1 


where the subscript zero indicates free-stream condi 
tions. The continuity equation is automatically satis 
fied through use of a stream function, and only Eqs 
(2.6) to be considered. These can 


(2.2) and remain 


be transformed to the (x, 7)-plane, using the trans 


formation formulas, 


[0,/ Ov], = (1/2) V potto/pox (0/0) ] 


[d dx], = (0/dx) — (1/2)x-'9(0/dn) \ 


When the given transformation formulas are applied 
to Eq. (2.9), there results the following expression for 


the vel CIty components: 


u = (1/2)uo(po, p) (df an 


a l atom) df r) (2. Ig 
oN pox \p / " be ; 
Using Eqs. (2.12) and (2.13) and dropping ter.ns that 
involve derivatives with respect to x at constant n, one 
can transform the x-momentum and diffusion equations 
into the following ordinary differential equations, re- 
spectively : 


f(g 'f')' + [bilge 'f')')’ = 0 (2.14 


AMN10" a QD» fino! — 0 (2? 1d 


where 


1 1 


od: = wo ‘uw; db = po 'p; A= vw Wry 


and primes denote differentiation with respect to 7. 
Apart from the strongly nonlinear character of Eq 
(2.14), it should be observed that the solutions are re- 
quired to satisfy two-point boundary conditions —1.e., 
conditions at 7 = 0 and » = ©. The method of solu- 
tion most suitable to problems of this kind is one that 
has been utilized in several papers dealing with bound 
ary-layer theory. Notable among these, perhaps, is 
the early paper by H. Schuh,* which describes the given 
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method of solution and also gives some consideration 
to the diffusion problem. The essential steps of the 
procedure are, first, to transform the differential equa- 
tions to a set of integral equations; and, second, to solve 
the latter by the method of successive approximations 
the specific advantages of which will be more fully 
discussed in Section (3) ]. 

Eqs. (2.14) and (2.15) are easily converted into the 
required integral forms. Eliminate the square brackets 
of Eq. (2.14) by carrying out the differentiation, and 
divide both sides by (q@:~'f’)’. 
take the antilogarithm of both sides. 


Integrate once and 
After perform- 
ing two further successive integrations, one obtains 
the following integral equation corresponding to Eq. 
2.14 


= [fh Ky foe fo x 
0 J0 
j 


exp ton / 1 fdn dyndn (2.17) 
/0 
The integration constant A, can be adjusted through the 


condition that lim f’ = 2, which means, according 


to Eq. (2.13), that the x-component of velocity ap 
The 


integration constant |f]o, which according to Eq. (2.13) 


proaches the free-stream value asymptotically. 
is related to the normal velocity component at n = 0, 


may be left unspecified. The third integration con 





20) 


16; iu 


8+ rd 
/j’— BLASIUS SOLUTION 
 ..-\st ITERATION 
—- 2nd ITERATION 
—- 3rd ITERATION 


—-- 4th ITERATION — 








6t 3rd & 4th 
Sy / ITERATION 


Ss 
> 
Sy 








0 + : + 

) rl «.« 3 

Convergence of successive iterations: 
hydrogen diffusion field 


isothermal case 


BOUNDARY 


f =f’ =0; mo 


LAYER 157 


stant has been set equal to zero, in accordance with 
the no-slip condition that is generally imposed upon the 
v-component of velocity at the wall. 

In a similar manner one can derive the following in 
tegral form of the equation of diffusion 


K wets 2k; 1 fdn bdy (2.18 
[miolo + Az exp )— od: 'fdn ¢dn (2.1 
70 od 0 


d) Sample Problems 


Nio = 


To illustrate the convergence of the successive ap 
proximations and to gain some insight into the nature 
of a diffusion field and its effect on the velocity profile, 
Eqs. (2.17) and (2.18) have been solved for two special 
cases: diffusion of hydrogen and helium into atmos 
spheric air at standard (N.T.P.) thermodynamic con 
All requisite physical data could be obtained 
with the ex 


ditions. 
International Critical Tables,‘ 
The latter was 


[see 
| See 


from the 
ception of Dj» for helium-air mixtures. 
estimated from some data given in reference 8 
Section (:3) |. 

According to Avogadro's law, the density of a binary 
gas mixture is a linear function of the composition pa- 
rameter. The same is not exactly true for viscosity 
coefficients and other transfer characteristics, but ex 
that the deviations from 
For the purpose of this 


perimental data indicate 
linearity are generally small. 
investigation it will be sufficiently accurate to represent 


the functions ¢,, in the form 


d, 1 + Bn, 2.19) 


where the B, are constants. 
From the appropriate experimental data, one can 
derive the following set of constants characterizing 


the two given problems: 


HI, — Air He — Air 
sj, = —0.507 0.097 
Bo = —0.931 —(). 862 
A = 1.63 3.94 


The subject of boundary conditions will be discussed in 
Section (4). Here we state only that the boundary 
conditions imposed in the given sample problems are 


as follows: 


laty =0; f’ =2; n 0 
at y (2.20) 

The zero-order solution for f, which may be taken to 
correspond to the case where foreign gas is entirely ab 
sent from the boundary layer, turns out to be the well 
known Blasius function. The zero-order solution for 
Ny, Corresponding to the same case, is identically zero. 
Convergence of the process of successive approxima 
tions, as illustrated in Fig. 1 for the hydrogen diffusion 
problem, is entirely satisfactory. The numerical cal 
culation can be greatly shortened, of course, by antic 
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ipating the course of future iterations in order to pick 
a better trial function. 

The velocity and concentration profiles for both hy 
drogen and helium diffusion are given in Fig. 2. It is 
interesting to note, in the first place, that the concen 
tration protiles are almost identical. This fact is indic 
ative of a certain insensitivity of the concentration pro- 
files, which is also confirmed by the rapid convergence 
indicated in Fig. |. So far as the velocity profiles are 
concerned, we can conclude (at least for the isothermal 
case considered here) that low density tends to lessen 
the velocity gradients, while low viscosity has the op 


posite effect. 
(3) THE GENERAL THEORY 


(a) Basic Equations 


The fundamental boundary-layer equations for vari- 
able fluid properties and for vanishing external pressure 


gradient may be given in the form* 


Ou Ou O Ou 
pu +o = (3.1 


Ov Oy 7 oy Oy 
Op ov = 0 (3.2) 
(O/Ox) (piu) + (0 Ov) (pv) 0 (33) 
O O O oT ou\? 
Uu (c,1) + pv (cpl) = (« ) ( ) 
x. wit ten ov\" av) * “Nay 
(3.4) 


The first two equations express conservation of .- and 


y-momentum, respectively. The third equation ex 
presses conservation of mass, while the last equation ex 
presses conservation of energy. The quantities p, u, 
k, and c, depend on temperature and on the concen 
tration parameter. For this reason it is necessary to 
introduce an additional differential equation describing 
the diffusion field before the given system can be inte 
grated. 

The diffusion equation that is applicable to the gen 
eral case contemplated in this section is not identical 
with the diffusion equation for isothermal flow, the 
latter being a special case of it. This circumstance 
arises, in the first place, from the dependence of the 
total particle density on temperature and, in the 
second place, from the fact that the diffusion current 


depends not only on concentration gradient but on tem 


perature gradient as well. Let C; represent the aver- 
age peculiar velocity of molecules of type 7, as defined 
by Chapman and Cowling.» We define the diffusion 


— 


current /, to be the number of foreign gas molecules 


jie 
flowing across a unit area perpendicular to C; in unit 
time—viz., 

J; = mC, ia.) 


* Neglecting diffusion thermo effect. 
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Effect of hydrogen and helium diffusion fields in air 
isothermal case 


Fic. 2 


where , represents the particle density of foreign 
molecules. Conservation of foreign molecules implies 
that the number of foreign molecules entering a given 
volume in unit time through convection equals the 
The 
form of the diffusion equation is, therefore, 


number departing through diffusion. general 


u(On,/Ox) + (OM OY) —div J, (3.6 


The variability of 7 with temperature, as mentioned 
above, implies that it is no longer possible to replace 
n, by the dimensionless parameter nj) = ,/n in the 
convection term. 

From the formulas given in reference 5, one cat 
derive the following expression for the diffusion current :7 

Jy. = —nDywd grad nw — ee al r\ (37 
‘. sa f 
where Dj» again represents the coefficient of mutual dif 
fusion and a@ is the thermal diffusion constant. The 
second term in Eq. (3.7) represents an additional dif 
fusion current caused by the temperature gradient and 
expresses the tendency of lighter molecules to migrate 
into regions of This 
known as “thermal diffusion,’’ and the thermal diffusion 


high temperature. process 1s 


constant is known for several binary gas mixtures.° 


Substituting Eq. (3.7) into (3.6) and neglecting de- 


rivatives with respect to x, as usual, one obtains the 


Tt Subscript 1 refers to the lighter molecules (cf. reference 6). 
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final form of the diffusion equation for the given prob- 


lem, 

on On, O | On, QN Hoof OL 
u Te sad nD» 7 ae ( 

On Oy oy | Oy 1 oy/ Js 


b) Method of Solution 


The problem at hand consists in finding a solution to 
the system of differential equations constituted by Eqs. 
3.01)—(3.4), 
given set of boundary conditions. 


together with Eq. (3.8), which satisfies a 
Eq. (3.2) has al 
ready been utilized through the implicit assumption 
that the pressure is everywhere constant and can there 
fore be neglected. The equation of continuity will be 
satisfied automatically through use of a stream function 
and can also be neglected. 

We begin, as in Section (2), by introducing a stream 
function y, a new independent variable 7, an f-function 
that is analogous to the Blasius function, and some di 
mensionless parameters through the following defini 


tions: 
py OW py OW ; 
ul : v= (3.9 
p OV p Ov 
l Poll - » 
n = \ : Y= V yoru f(n (3.10 
Zz Mod 
Mm imo; p = depo; k = gsk 
Cy = Opi NDi2 = dso; 1 den (3.49) 
T= tT 


Again we transform the differential equations to the 
v, n)-plane by means of the transformation formulas. 


2 | — 1 g/ porto 0 | 
ov |, oW ae on 


| Oo | (a) | O 
= — Y 'n 
Ox_], Ox 2 On 


One finds, as before, that the velocity components are 
given by the expressions 


u = (1/2)uo[go.—'f"] | 


1 . | uotto 


a 


Using this result, together with the transformation 


-<" ‘iy’ - Al) 


> 


formulas given in Eq. (3.12), and setting partial de- 
rivatives with respect to x at constant » equal to zero, 
one obtains the following system of ordinary differential 
equations, expressing conservation of x-momentum, 
conservation of energy, and conservation of foreign 
molecules, respectively : 


wd dg 


Sloe 'f')' + [b1(¢2-1f) 


flout)’ + Pro 1 et’ |’ + ] 4( vo — 1) x 


M7: [(¢2'f’)’}? = O (3.15) 
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QD» Mf [dem y + [b5( 1) ‘ QaN;(1 — Nyt 1’ |’ 0 


Prandtl Mach Numbers 


have been introduced through the identities 


where the free-stream and 


Pro = Cymo/ R / 
) 17 
(44 1) My? = mo? Acp,To\ 
The next step will be to transform Eqs. (3.14 


3.16) into integral equations that lend themselves 
readily to the process of successive approximation 
Variation of 


This is accomplished by the ‘“‘method of 


parameters ‘as described, for instance, by Ince.‘ Let 


there be given a second order nonhomogeneous linear 


differential equation of the form 
Lin\u } 


and let the solution to the associated homogeneous 


u U,lo t A / WV, dn 


equation be 


and A, are constants of integration. It 


where [z, | 
can be shown, in this particular case, that the solution 
to the complete equation has the form 


ul | 24 lo + AK / WV dy T 
/ V | ae in| dn 53.1S 


This mathematical theorem is readily applicable to the 


given problem. Since the resulting system of integral 
equations will ultimately be solved by the method of 
successive approximations, there is some freedom in 
selecting terms to be included in the general term F 
representing the nonhomogeneous component of the 
differential equation. As will be seen below, the F 
terms appearing in the present development include 
mainly terms that are peculiar to the diffusion problem, 
while the W-functions correspond essentially to the 
terms considered in standard theories of the compres 
sible boundary layer. 

Let the variable uw, occurring in Eq. (3.18) represent 
1, 2, and 3, 


the quantities f, /, and my) for n = respec 


tively. By means of elementary calculations one ob 
tains the following expressions for F,, and W,,: 


v, = 6. | o ' exp ) -{f ?) if dn sdn| 
Jo ( So (3.19) 


F, = 0 | 
oe ae hee 
WV. = ¢$;~! exp y—f rn 3 | dsf dn 


F, = — Props} (1/4) (vo — 1) 
Mooi [(¢2~'f")’]? + fos’t} 


(3.20) 
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J0 . 9 
F; ~h5 |} oo 'he’ fry 


lads, (1 


The given set of integral equations is amenable to 
solution by the method of successive approximations. 
One of the advantages of the present formulation lies 
in the fact that explicit constants of integration are 
available which can be readily adjusted to satisfy the 
two-point boundary Other 
are the comparative ease with which the requisite nu- 


conditions. advantages 


merical operations can be carried out and the possibility 
The 


calculation can be performed by a subsidiary worker 


of detecting and controlling numerical errors. 


with the aid of a manual calculating machine and re- 


quires approximately 10 hours for each iteration. 


c) Sample Problems 


The sample problems considered in this section were 
chosen to illustrate the effect of a diffusion field of helium 
atoms on the dynamic and thermal characteristics of a 
compressible boundary layer at high Mach Numbers. 
To render the solutions fully determinate, the following 
conditions must be specified: (a) boundary condi 
tions; (b) exact functional dependence of ¢, on f and 70, 
and the numerical values of P7o, Jo, yo, and a. 

The following boundary conditions are assumed: 


f=f' =0; t=6; m=1 aty=0 , 2 


f =2; t=1; mo=0 aty= a 


In particular, it is supposed that the concentration of 
helium near the wall is constant, while the normal com 
ponent of velocity vanishes there. Possible objections 
to this supposition, based on the question as to whether 
conditions of this nature are physically feasible, will be 
answered in the next section, where the practical sig 
nificance of the given theory is discussed. 

To facilitate the numerical calculation, the @-fune 
tions are represented as a product of two factors, de 


vending, respectively, on 7) and ¢, 
I 2 
PD» = Dri( No) Py2(t) (3.2 


The functions ¢,:, which express dependence of a given 
property on the local concentration of helium, can be 
represented with sufficient accuracy by linear expres 
sions. The constants entering into these linear ex 
pressions are assumed to be independent of temperature 
The 
resulting functional relationships are given in the fol 
lowing list, together with the sources from which rele- 


and are usually evaluated at N.T.P. conditions. 


vant empirical data were obtained: 


on l O.09TH Reference 4 
1.505 ; 
ty re Sutherland law 
1 + O.505t7! 
de I] O.862n Reference 4 


o s {= Gas law 
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Fic. 3. Effects of a helium diffusion field: 17) = 12; [tf 6 

rei) | 1.S7n Reference 4 

30 = 1 Assuming constant ¢ ind con 
stant Prandtl] Number 

ri) = ] 4.18 Reference 5 

oy = 1 Usual approximation 

o l n and Dj. both independent of 
concentration 

o 2.98¢ Estimated from data for // A 
He-A and Hy:-air mixtures giver 
in references 4 and 8, and as 
suming a free-stream  temper- 
ature 7) = 392°R 

ou = | Avogadro's law 

de @ Gas law 

We further assume 

Pr 0.75, vo = 1.40, a = 0.365 


Only the free-stream Mach Number remains now t 


be specified. Two sample calculations have been 
carried out, corresponding to the conditions J/, 12 
and S, respectively. The resulting velocity, tempera 
ture, and concentration profiles are indicated in Figs. 
5 and 4, together with the corresponding velocity and 
temperature profiles that obtain in the absence of a 
diffusion field and which have been calculated by Van 
Driest.? The theoretical and practical significance of 
these results will be discussed in the following section 
Here we shall only state that approximately three iter 


ations are required to achieve sufficient convergence, 
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provided one uses Van Driest’s velocity and temper 
iture profiles for the first trial functions. Faster con 
vergence can of course be achieved through the use of 


more suitable trial functions 


{) RESULTS AND APPLICATIONS 


a Interpretation of Results 


The most striking effect of the helium diffusion fields 
is observed in the temperature profiles, as indicated by 
Figs. 8 and 4. The normal peaks in the temperature 
profile, which obtain whenever the wall temperature 
is kept at a value less than adiabatic or equilibrium 
temperature, are caused by the dissipation of kinetic 
energy. It is evident that these peaks are greatly 
diminished by the admixture of helium, and in par- 
ticular instances they may be altogether eliminated. 
[his result is mainly due to three contributing factors, 
the first two ot which act to reduce the temperature rise 
caused by dissipation of kinetic energy, while the third 
factor tends to reduce all temperature gradients in the 
boundary laver. It was pointed out in Section (1 
that helium is notable for its low density, high specific 
heat, and high heat conductivity. The first of these 
properties tends to reduce the velocity gradient in the 
boundary laver, as was shown by the results obtained 
in Section (2). In view of the fact that the dissipation 
of kinetic energy varies as the square of this gradient, 
one can conclude that low density tends to reduce vis 


cous heating. For a given rate of dissipation of ki 
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Fic. 4. Effects of a helium diffusion field: J/) = 
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netic energy, the resulting temperature rise will be re 
duced by the presence of helium, because the latter has 
Finally, all 


gradients in the boundary layer tend to be reduced on 


a very high specific heat. temperature 
account of the high thermal conductivity contributed 
by the helium atoms. In this manner each of the three 
most unusual physical properties of helium gas plays 
a part in reducing the peaks and gradients of hyper 
sonic temperature profiles. 

The two sample problems of Section (3 
While the magnitude of the initial 


differ in the 
following respect 
gradient of the temperature profile at 1/ 12 is 
greatly reduced by the admixture of helium, its sign 1s 
unchanged. This means that heat still flows into the 
wall, though the rate at which this takes place may be 
greatly reduced. At .l/ S, however, the sign of the 
initial temperature gradient has changed, causing a 
reversal in the direction of heat flow. These results 
are easily interpreted by means of the well-known ex 
pression for adiabatic wall temperature, which is ap 
homogeneous boundary 


plicable to any chemically 


layer with Prandtl Number approximately equal to 


To = Tofl + (12) (4 l)a/Pr M 1.1 


This formula enables one to express the adiabatic wall 
temperature of an all-helium boundary layer in terms 
of the adiabatic wall temperature in air, corresponding 
to the same free-stream temperature and velocity 
For this condition (4 1).17° varies inversely as thy 


specific heat, which implies, according to Eq. (4.1), that 


l To) ye 
1 To) air 


(VY Pre, He Cp) ai 
: eS 0.193 (4.2 


V Pre, siz Cp) He 


Applying this result to the first problem, for which .1/ 


= 12, one finds that the adiabatic wall temperature 
would be reduced from its value 25.97% in air to 5.827 


Assuming an “‘aver 


in an all-helium boundary layer. 
age concentration’’ of helium equal to SU per cent, the 
estimated adiabatic wall temperature becomes 7.027 

which is greater than wall temperature for the given 
problem. This explains the small but positive tem 
perature gradient at the wall which has been calculated 
At a free-stream Mach Number of 8, 
adiabatic wall temperature is reduced from 12.17) in air 
to 3.147) in an all-helium boundary layer. Assuming, 


however, the 


as before, an average concentration of helium equal 


, which 


This 


to SO per cent, this value is increased to 3.687 
is stall considerably less than wall temperature 
explains the reversal in heat flow noted above. 

One further point of interest concerning the temper 
ature profiles must be explained. It is somewhat 
remarkable that there could occur a pronounced dip 
in the temperature curve near the outer edge of the 
boundary layer, where temperatures below the free 
stream value are maintained. At least to the present 


author’s knowledge, such a phenomenon is not indi 
cated by other theories of the laminar boundary laver, 


but its appearance in the present context is readily 
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plainable. An inspection of the numerical calculation 
discloses that the phenomenon is caused by the term 
(fes’t) occurring in Eq. (3.20). It arises from the con- 
vection term of the original energy equation and ex- 
presses the effect of variable specific heat on convection 
of enthalpy. From the given concentration profiles 
it is evident that strong negative gradients of specific 
heat exist near the temperature minimum. The fol- 
lowing simplified argument indicates how these gra- 
dients account for the unusual behavior of the tem 
perature profiles: Consider a small rectangle fixed in 
the outermost region of the boundary layer, where the 
velocity gradient is effectively zero. Two of its sides 
are taken to coincide with the local stream lines, and 
the remaining sides are perpendicular to the flow direc- 
tion. Let the upstream and downstream sides be re- 


ferred to as A and B, respectively. A fluid particle 
traversing the given rectangle moves in the direction of 
Suppose the temperature in the rectangle 


Then, be- 


decreasing 7. 
were maintained at its free-stream value. 
cause of the negative gradient of specific heat, the fluid 
particle would gain enthalpy while traversing the rec- 
tangle. Since dissipation of kinetic energy can be 
neglected in this region, the temperature within this 
rectangle will naturally drop. In this way a positive 
temperature gradient is set up which opposes the con- 
vection of enthalpy caused by the variation of specific 
heat. 

Two remarks should be made in conclusion. A com- 
parison of numerical values indicates that the thermal 
diffusion term contributes only a small correction to 
the concentration and could be neglected in the given 
instances without causing undue errors in the result. 
This is partially due to the fact that the presence of 
helium in the boundary layer has greatly reduced the 
strength of temperature gradients, thereby causing the 
thermal diffusion term to be smaller than an initial 
estimate would disclose. Further, it appears from the 
figures that the velocity profiles are not greatly affected 
by the admixture of helium, although the initial gra- 


dients are somewhat reduced. 


b) Concerning the Boundary Conditions 


A few words should be said concerning the particular 
boundary conditions that have been assumed in the 
sample problems. The objection may be raised that it 
is not physically feasible to maintain a constant non 
vanishing concentration of foreign gas near the wall 
without maintaining, at the same time, a nonvanishing 
normal velocity through the wall. 

One may reply, in the first place, that questions re- 
lating to the physical feasibility of the given boundary 
conditions need not concern us at all. They have been 
selected for the purpose of investigating the effect that 
the presence of a foreign gas can have on boundary- 
layer profiles, unmasked by any additional disturb- 
ances, such as the existence of a nonvanishing normal 
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velocity near the wall. If this condition should be a; 
idealization that is physically unattainable, it never 
theless describes trends and tendencies of general sig 
nificance. Moreover, the assumption of vanishing 
normal velocity near the wall is accidental not only t 
the effect of diffusion fields on the boundary laver but 
also to the theory and method of solution presented in 
this paper, which is directly applicable to the more gen 
eral case. 

On the other hand, the conditions that have beer 
supulated in the sample problems are by no means ex 
perimentally inconceivable. Let us suppose the wall 
to be a surface that is perfectly permeable to both pri 
mary and special molecules. Such a surface would pro 
hibit the existence of a discontinuity in composition 
If the gas pressure inside the surface is maintained at 
the free-stream value, there will be no net flow of gas 
through the wall. In this case the concentration 
special molecules in the boundary layer 1s maintained 
by means of a pure diffusion current through the 
wall, in keeping with the stipulated boundary condi- 


tions. 


c) Applications and Conclusions 


The theory and results presented in this paper have 
bearing on the problem of control and reduction of aero 
dynamic heating at hypersonic speeds. The special 
feature of the present approach to this problem lies 
in the utilization of diffusion fields for the purpose of 
reducing the detrimental effects of viscous dissipation, 
a procedure that may be compared with the applica 
tion of a lubricant between two surfaces. Although the 
given theory and method of solution are adapted to a 
fuller investigation of this problem, the numerical ex 
amples presented in this report are mainly concerned 
with diffusion fields of helium, with which highly satis 
factory results have been achieved at the free-stream 
Mach Numbers of S and 12. Whereas at the higher 
Mach Number the influx of heat was practically elim- 
inated, a reversal in the direction of heat flow has been 
effected at the lower Mach Number. 
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Problems and Results of Investigations on 


Cascade Flow ™ 


HERMANN SCHLICHTINGT 


Institute of Flutd Mechanics, Technical Unwersity of Braunschweig, Germany 


SUMMAR\ 


A summarizing report is given of the research work on the 
flow through cascades of airfoils which has been conducted in 
recent Vears at the Institute of Fluid Mechanics at the Technical 
University of Braunschweig. The theoretical investigations led 
to a simple solution of the ‘Indirect Problem” and the ‘Direct 
Problem” for the two-dimensional incompressible frictionless 
flow through cascades. By an extensive computation of ‘‘Cas 
cade Downwash Tables,”’ the amount of calculation for any spe 
cial case has been reduced to such an extent that a systematic 
theoretical investigation of the cascade problem has become pos 
sible 

Two-dimensional viscous flow through cascades has been 
treated by applying boundary-layer theory, thus giving for the 
first time theoretical data for the loss coefficients. An extensive 
theoretical program, which aims to explore the relations between 
the geometric and aerodynamic parameters of cascades, including 
loss coefficients, has been started and partly completed 

For the two-dimensional compressible flow through cascades at 
high subsonic velocities, a convenient method of computation 
has been given which is based on the Prandtl-Glauert rule In 
this way, the above-mentioned solution of the ‘Direct Problem”’ 
for incompressible flow has become applicable also to compressible 
subsonic flow 

The theoretical results have been checked by some systematic 
experiments on two-dimensional cascades. Special attention has 
been given to obtain two-dimensional flow by means of boundary- 
layer slots at the walls and by applying suction. Pressure dis- 
tributions, as well as loss coefficients, as obtained from experi- 
ments are in good agreement with theory 

Furthermore, numerous experiments have been done with 
cylindrical axial-flow cascades having various hub ratios and 
untwisted blades These cascades were tested with axial inflow, 


1 considerable swirl of the outflow. By comparison 


thus giving 
with tests on two-dimensional cascades with blades of the same 
profile, it has become evident that for cylindrical cascades the 
effects of radial divergence—i.e., the variation of the solidity 
ratio over the radius—on the pressure distribution over the blade 


ind on the local loss coefficients generally are extremely small 


Received July 1, 1953 Revised and received September 24, 
1953 

* From a lecture held at the Conference of the Committee for 
Fluid Motion Research of the ‘‘Verein Deutscher Ingenieure”’ 
in Karlsruhe, Germany, on March 19-21, 1953. The full German 
text has been published in Zeitschrift fiir Flugwissenschaften, Vol 
I, No. 5, pp. 109-122, October, 1953. For this English trans 
lation I am indebted to Dipl.-Ing. E. S. Krauss, AFRAeS, Lon 
don, and to Prof. J. Weske, of the University of Maryland 

+ These investigations have been conducted with the assistance 
Forschungsgemeinschaft,’’ Bad Godesberg, 
and of several German industrial firms. Among the staff of the 
Institute, Dr.-Img. N. Scholz, Dr.-Ing. L. Speidel, Dipl.-Ing. H 
Schaffer, and Ing. H. Goldmann prominently participated in the 
work 


t Professor. 


of the ‘‘Deutsche 


Also, some research has been done on two-dimensional cascades 


with clearance 


(1) INTRODUCTION 


F emia UL INVESTIGATION of the cascade prob 
lem only by systematic experimental research, 
as it was done for the single airfoil by extensive wind 
tunnel tests in U.S.A. during the last decades, seems 
to be quite hopeless in view of the much greater number 
of parameters and also because of heavy experimental 
difficulties. Real progress beyond occasionally meas 
ured special cases can therefore only be achieved by 
investigations, which, however, 


extensive theoretical 


require a careful supervision by experiment. From 
these considerations the Institute of Fluid Mechanics 
at the Technical University of Braunschweig decided 
to carry out closely interrelated theoretical and ex 
perimental investigations. 

The leading idea of our research work is that a deepe r 
knowledge of the flow phenomena, as they occur in 
turbomachines in a most complex manner, will only be 
possible by an extensive analysis. This analysis must 
extend to theory, as well as to experiment. The basis 
of all investigations of the flow through cascades is the 
incompressible flow through a two-dimensional cas 
cade. Only if this basic problem is fully understood, 
will it be profitable to attack the problems involved in 
the remaining parameters. 

The most important of these other parameters are 

(1) Three-dimensional effects (a) along the walls and 
(b) due to clearance. 

(2) Radial divergence of cylindrical cascades.** 

3) Compressibility effects (Mach Number 

(4) Centrifugal forces (for the rotating unit 

Systematic research on incompressible flow through 
a cascade seems to be unavoidable if a sound knowledge 
of the physical process is to be acquired instead of being 
satisfied with purely empirical results. 

Therefore, at the beginning, we have confined our in 
vestigations mainly to the incompressible flow through 
two-dimensional cascades, since so many unsolved 
problems exist in this particular field. 

Among the other parameters mentioned before which 


complicate the problem of cascades, up to now we have 


** Variation of the solidity ratio with radius. 
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Geometric and aerodynamic parameters of a two 
dimensional cascade. 


Fic. 1 


only dealt with a few and to a moderate extent. There- 
fore I should like to report at first on the investigations 
of two-dimensional cascades and later on some of the 
in particular, the effect of Mach 
especially end 


other problems 


Number; three-dimensional effects, 
losses and radial divergence as occurring in cylindrical 


cascades. 


(2) INVESTIGATIONS OF TWo-DIMENSIONAL CASCADES 


Whereas in earlier papers the flow through cascades 
has been treated as a flow through a passage, all mod- 
ern investigations are based on the concepts of airfoil 
theory. Each blade is treated as an airfoil and the 
mutual interference of the blades plays a most im- 
portant role. The main problem in investigations of 
the two-dimensional cascade is to find the relation be- 
tween the geometric and aerodynamic parameters of the 
geometric parameters are, as 


cascade. The main 


shown in Fig. |, the profile of the blade: 


with the maximum camber, fc 
thickness chord ratio, ¢ ¢ 
camber line, y,(v) 

thickness distribution, y,(v) 


and the geometry of the cascade as given by: 


solidity ratio, dc 

angle of stagger (or blade angle), 8 
The first parameters that define the profile of the blade 
are known from the airfoil alone. The two additional 
parameters defining the geometry of the cascade depend 
on the arrangement of the blades in the cascade. 

The aerodynamic parameters of the cascade are: 


inflow and outflow velocity, wy, w» 


ICAL SCIENCES MARCH, 1954 
angle of inflow and angle of outflow, 3,, 32 or aj, a 
deflection of flow deviation, AB = 38, — 23 
pressure difference across the cascade, 1p = p» bp 
chordwise pressure distribution over the blade, 
resultant aerodynamic force on the blade, R( LU, S 
with the tangential component, / 
and the thrust component, S 


loss of total head in the cascade, A/ 


The establishment of the relation between these 
geometric and aerodynamic parameters is the aim of 
these investigations. This relation is greatly affected 
by the Reynolds Number and Mach Number. For 
the single airfoil, the relation between the geometric 
and aerodynamic parameters has been established, to 
a large extent, by aeronautical aerodynamics. But 
with regard to cascades, our knowledge is far from being 
satisfactory. The reason is that, because of the greater 
number of geometric and aerodynamic parameters, the 
problem of the cascade is much more difficult than that 


of the single airfoil. 


2.1) Theoretical Investigations 


In all previous theoretical investigations of the flow 
through cascades a frictionless fluid has been assumed 
throughout.' From the very beginning of our investi- 
gations, we have taken into account the friction of the 
fluid for the first time.” * Thus, a theoretical calcu- 
lation of the loss coefficients has become possible. This 
is of particular importance because these losses deter- 
mine the quality of the cascade and thus the efficiency 
of the turbomachine. 

(2.1.1) Incompressible Flow.—But before we could 
start with the investigations of the cascade flow with 
friction we had to deal with extensive basic problems 
in the field of frictionless flow through cascades. 

For the calculation of the boundary laver, it is neces- 
sary to know the pressure distribution according to 
potential flow theory. Unfortunately, the existing 
methods available for calculating this pressure dis- 
tribution proved to be unsuitable for our purpose of 
systematic investigations of the cascade flow, since 
they are either too elaborate or they do not satisfy one 
requirement which is to be able to vary the geometric 
parameters of the cascade separately. The possibility 
of varying the geometric parameters separately —e.g., 
variation of solidity ratio with the remaining param- 
eters—the blade angle and the blade profile being con- 
stant, proves to be of predominating importance in 
order to obtain clear information on the relation be 
tween the geometric and aerodynamic parameters. 

The two most important basic problems of the theory 
of cascades in frictionless fluid are the so-called ‘‘In- 
direct Problem” and the “Direct Problem.” 

In the ‘‘Jndirect Problem’ the known quantities are 
the directions of the inflow and outflow—i.e., the 


velocity diagram—and the unknown quantities are the 
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geometry of the cascade and the geometry of the blade 
section and also the pressure distribution on the blade. 

In the “Direct Problem’ the known quantities are 
the geometry of the cascade and the geometry of the 
blade profile, and the unknown quantities are the re- 
sulting aerodynamic force on the blade, the pressure 
distribution and the direction of the outflow—all of 
these as functions of the direction of the inflow. 

In the ‘Indirect Problem" the aerodynamics of the 
cascade are given and the geometry is asked for. 
“Direct Problem"’ is the reverse of the ‘“‘Indirect Prob- 
lem.” 

We have succeeded in finding a convenient method 
for solving these two problems. The solution is based 
on the method of singularities, which replaces the blade 
distribution of sources, sinks, and vortices—a 
method already suggested earlier by Schilhansl* and 
Ackeret.’ A solution for the “Indirect Problem” has 
been given by Scholz'’ and a solution for the ‘Direct 


by a 


Problem” by Schlichting. '® 

With the aid of a catalog that compiles universally 
tabulated functions, the amount of computational work 
has been reduced considerably as compared with the 
work involved previously. For instance, the work re- 


quired for the “Direct Problem” of a given cascade 


The. 


with ten different angles of inflow can now be com- 
pleted in about 20 hours instead of 20 days.* * Hence, 
a systematic theoretical investigation of the cascade 
problem seems to have become possible and we have 
started with it. 

A primary result of such a systematic computation 1s 
shown in Fig. 2. The lift coefficient c, is plotted against 
the angle of inflow a, for a series of cascades with 
blades of the NACA 0010 profile. 
are dc = 0.5 and 1.0 and the blade angles 8g = 90°, 
60°, 45°, and 30°. 
single airfoil, dc = @, 
the well-known fact that the slope of the lift curve of 
the blade in cascade is only a fraction of the lift slope 


The solidity ratios 


The corresponding values for the 
are included, thus confirming 


of the single airfoil. 

It is of particular interest to note that, except for the 
case of the unstaggered cascade, 8, = 90°, the curves 
of c, plotted against a; do not pass through the origin 


This means that, at an inflow parallel to the chord, 
a, = O, the lift coefficient is different from zero. Ina 
staggered cascade of symmetrical profiles of finite 


thickness with inflow parallel to the blade chord, the 
flow is deflected toward the cascade axis, which was 
found first by Ruden* and later by Marcinowski.’ In 
Fig. 2c, the angle of deflection a2 is plotted against the 
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Fic. 3. Pressure distribution of turbine blade sections with large camber in two-dimensional cascades; theory according to Richter ;™ F 
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Fic. 4. Unstaggered cascades: 8g = 90°. Thickness distribution: NACA standard four digit airfoils, camber line—parabolic. 





fc TI 
wit 





*hter :* 


bolic. 





INVESTIGATIONS 















































FIG 5 
of maximum camber f/c, according to Speidel.!” 
boundary layer fully turbulent. @C 
of entry 


blade angle 3, for various values of the solidity ratio 
d c for an inflow parallel to the chord of the blade. 

In the course of our systematic investigations on cas- 
cades we have made extensive calculations of pressure 
distributions, because these are the basis for all theo- 
retical computations of the coefficients. Our 
method of calculation is most suitable for blades of 
As, how- 


loss 


moderate thickness and of moderate camber. 
ever, blades with large values of maximum camber 
occur in steam turbines, we have extended our method 
to such cases. In Fig. 3, the theoretical pressure dis- 
tributions of a turbine cascade, having blades with large 
values of maximum camber and large deflections of 
flow, are compared with experimental pressure dis- 
tributions for two different values of the solidity ratio. 
The agreement is good, particularly in view of the dis- 
crepancy one would have expected to exist between 
theoretical and experimental values due to the neglect 
of the effect of friction. The rapid change of curvature 
of the contour of the upper surface of the airfoil section 
considered is reflected in the pressure distribution in 
form of a sudden change of the pressure for the cascade 
Further 


with the smaller solidity. pressure distri- 





Theoretical loss coefficients ¢ according to Eq. (1) for unstaggered two-dimensional cascades of Fig 


Blade angle 6g 
first occurrence of separation on upper and lower surface, respectively d = ideal angle 
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4 for various values 
90°; thickness t/c = 0.15; Reynolds Number Re = wac/v = 10°, 


butions, together with their correlation with theory, 
will be given later. 

Viscous flow, loss coefficients: The main objective 
of our investigations is to find the most favorable cas 
cades with regard to flow energy. The 
method of calculating the losses is basically that of the 


losses of 


boundary-layer theory applied to the cascade. We 
have embarked on an extensive program on these lines, 
which is partly completed. 

In Fig. 4, a survey is given of a number of cascade 
configurations, the losses of which have been calculated 
by Speidel.'’ This first part of our program refers to 
unstaggered cascades only—1i.e., the solidity ratio being 
the only parameter of the geometry of the cascade. 
All blade sections are NACA sections, with varying 
thickness chord ratios and varying maximum camber. 
From the numerous results, only one diagram about 
loss coefficients may be given here, Fig. 5. One of the 
most important aerodynamic parameters of a cascade 
is the velocity difference Aw, in the tangential direc- 
tion across the cascade. It is desirable to design cas 
cades with the greatest possible deflection of the flow 


at a minimum of losses. In Fig. 5 the loss coefficient 
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£ = Ah/(p/2)weia? = Ah/Qoia (1) 


based on the theoretical outflow velocity weg is plotted 
against the velocity ratio Aw,/w,, which corresponds 
substantially to the angle of deflection of the flow Ag. 
Each of the three diagrams refers to a constant value of 
maximum camber and three different solidities of the 
cascade, d/c = 0.5, 0.75, and 1.00. 

For moderate deflections of flow, the minimum losses 
are associated with the largest values of cascade spacing, 
because the losses are essentially proportional to the 
number of blades for small blade loadings. For larger 
angles of deflection, however, a closer spacing is more 
favorable. The reason for this behavior is that, at 
greater angles of deflection, the lift loading of the 
individual blade is becoming too large. For excessive 
loadings, the flow separates and the losses increase 
rapidly. The curve referred to as ‘“‘Optimum’”’ is the 
envelope of the three dc curves and indicates the mini- 
mum of losses associated with a given angle of deflec- 
tion and solidity ratio. The symbols @ and O indi- 
cate the beginning of separation on the upper and lower 
surface of the blade, respectively. Only in the region 
between © and @ is the flow attached on both sur- 
faces. It is of particular interest to note that the opti- 
mum values (envelope) are throughout on the right- 
hand side of @, which means that the flow on the upper 
surface is slightly separated. This is in good agree- 
ment with experiments but is completely different 
from the behavior of the single airfoil. Whereas in 
the case of the single airfoil an increase in lift coeffi- 
cient is always associated with a considerable increase 
in drag particularly in the range of partly separated 
flow, in the case of the cascade it is possible to obtain 
a reduction of the losses by increasing the spacing which 
corresponds to an increase in blade loading at constant 
angle of deflection. The reason is that the increase in 
drag of the individual blade is more than compen- 
sated by the decrease in the number of blades. 

(2.1.2) Compressible Flow (Subsonic).—In the ma- 
jority of practical cases, the velocity of flow through 
cascades is so high that the compressibility effects 
cannot be neglected. The methods known for calcu- 
lating the frictionless compressible flow through cas- 
cades**~* are too complicated to be suitable for the sys- 
tematic investigations that are most urgently needed. 
It is therefore feasible also to apply the Prandtl-Glauert 
rule, which for the single airfoil has proved most satis- 
factory, to the cascade in the range of subsonic ve- 
The important assumption for the applicabil- 

that the deviations 
also remains valid for 


locities. 
ity of the Prandtl-Glauert rule 
from a parallel flow are small 

the flow through cascades. This assumption therefore 
corresponds to a restriction to small angles of deflection 
of the flow. For the present analysis, the most con- 
venient form of the Prandtl-Glauert rule is that for 
which the blade section is the same in the compressible 
flow and in the corresponding incompressible flow. In 
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this case, therefore, the reduction in density due to the 
expansion of the gas corresponds to a smaller spacing, 
According to 
flow 


which in turn causes greater velocities. 
Scholz,'* the the compressible 
through a cascade with a given spacing and a given 
stagger at a given Mach Number can be reduced to the 


calculation of 


calculation of an incompressible flow through a cas- 
cade with a different spacing and a different angle of 
stagger. For the unstaggered cascade, the relations 
are extremely simple. 

In this case a cascade in compressible flow at a 
Mach Number J = w..a.. and a solidity d, c corre- 
sponds to a cascade in incompressible flow with a 
smaller solidity ratio d, c, which is given by the equa 


tion 
d,/c = (d,/c) V1 — M? (2) 


The equations for the transformation of the velocities 


and pressures are 


d, d 

Ax a) = («, ) VYi— 7 3) 
( c / 
d ) d 

pis: =, My} Di <; )/ V1l-—- Mw (4) 


By virtue of reducing the compressible subsonic flow 
through cascades to an incompressible flow, the method 
for solving the ‘Direct Problem” in incompressible 
flow can be readily used. 

As an example, the velocity distributions of unstag 
gered cascades with NACA 0010 blade section are shown 
in Fig. 6, the inflow being parallel to the blade chord 
and therefore with zero deflection. 

This diagram permits direct information regarding 
the chordwise position of the point at which the velocity 
of sound is first reached with increasing free-stream 
Mach Number. The local Mach Number J/jo.4; = 
w a.., which is obtained by 

Ce a ee (5) 

a Ma Ba W a 

is shown for the single airfoil and also for the cascades 
with the solidity ratios dc = 1 and 0.75. Whereas 
for the single airfoil the velocity of sound is at first at 
tained at a Mach Number of J = 0.8, for the solidity 
of d/c = 1 the velocity of sound is first reached at a 
Mach Number of J = 0.75 and for the solidity of 0.75 
at a Mach Number of JJ = 0.72. The curves shown 
in Fig. 6 are only valid for completely subsonic flow. 
As soon as the velocity of sound is reached locally at 
Mach Numbers near ./ = | in the smallest cross sec 
tion of the passage between the blades, compression 
shocks occur, and the passage is chocked and super 
sonic velocities will occur aft of the smallest cross sec 
tion. In Fig. 6, the curves in the supersonic region are 
therefore shown in dotted These 
have been investigated more closely in theory and ex 


lines. conditions 
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Chordwise velocity distributions over unstaggered cascades of various solidity ratios d/c for symmetrical inflow. 
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Comparison 


of theory and experiment: O@ measured velocities along the upper surface and lower surface, respectively 


periment by Ackeret and Rott.'' It can be generally 
stated that a blade in cascades is more sensitive to 


compressibility effects than the single airfoil. 


(2.2) Experimental Investigations 


I shall now deal with our experiments on two-dimen- 
sional cascades. Apart from the extreme difficulties 
of experimental investigations of two-dimensional cas- 
cades associated with obtaining fully two-dimensional 
flow conditions, there is the additional fundamental 
difficulty of the immense amount of experimental work 
required if a clear concept of the effects of various pa- 
rameters upon the cascade is desired. 

Experimental investigations of two-dimensional cas- 
cades therefore should be carried out, in my opinion, 
first and foremost for checking systematic theoretical 
investigations apart from supplying practical data in 
special cases (e.g., roughness and special blade sec- 
tions). 

(2.2.1) The 
Our two-dimensional cascade tunnel, which was con- 


Two-Dimensional Cascade Tunnel. 
structed after some years of preliminary research by 
Dr. N. Scholz and Ing. H. Goldmann, is particularly 


suited for such experiments. The experience gained 


with other tunnels has, of course, been used, in par- 
ticular the experience—insofar as it has been pub- 
lished—of Erwin and Emery" in the NACA Labora- 
tory at Langley Field. 

In Fig. 7, a schematic presentation of this cascade 
tunnel is shown. The tunnel is of rectangular cross 
section and can be attached to the nozzle of a wind 
tunnel with a diameter of 1.3m. The main dimensions 
are shown in Fig. 7. The number of blades is seven; 
thus, according to our experience, the flow correspond- 
ing to an infinite number of blades is obtained for the 
middle blade. The following parameters can be set 
at any required value within a satisfactory range- 

solidity ratio, d/c 

blade angle, 6, 

inflow velocity, w; 

direction of inflow, 6; 

The following quantities are measured as functions of 
these parameters: 

(1) Pressure distribution on the blade (pressure 
taps), p(x). 

(2) Direction of outflow behind the cascade, 8». 

(3) Outflow velocity, we. 

(4) The total head behind the cascade, /io. 
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section NACA 0010; Re = wec/vy = 5X 105, fully turbulent boundary layer. 
lines represent theory by Speidel!®; and the circles with projected lines mark the position of separation point at x/c 


The distributions of the last three quantities have 
been measured in a plane behind the cascade parallel 
to the cascade axis extending at least over one spacing. 
A special instrument was developed which made it 
possible to take the readings of all quantities in one 
survey, thus reducing considerably the time required 
for the tests. The greatest inflow velocity of 60 m per 
sec. is in the range of incompressible flow. The chord 
of the blade is c = 200 mm., and the Reynolds Number 
based on the outflow velocity is 


Re = we-c/v = 5 X 10° 


Special care has been taken to obtain two-dimensional 
flow through the cascade, which is rather difficult be- 
But the 
verification of a good two-dimensional flow is of ex- 


cause of the boundary layers on the walls. 


treme importance, particularly with respect to the 
comparison of theory and experiment. 

In order to remove the boundary layer, slots are 
arranged at the four walls ahead of the cascade which 
are used for removing the boundary layer by suction. 
The cascade tunnel with this arrangement for suction 
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Solid lines represent experiments by Scholz; dashed 
= 0.95 


is shown in Fig. 8. The two slots of each sidewall are 
connected by means of two ducts and each slot of the 
upper and lower wall by means of one duct, a total of 
six ducts in form of flexible tubes with an axial blower 
With 
this two-dimensional cascade tunnel we have obtained 
a large number of results, some of which I will now dis 


with an installed power of approximately 33 hp. 


cuss briefly. 

(2.2.2) Some Experimental Results as Compared with 
Theory.—The first series of tests was done with cas 
cades of blades of NACA 0020 profile. Though the 
thickness chord ratio of this section is comparatively 
large, (/c = 0.2, making this profile of limited practical 
value, it was chosen in order to test the interference 
effects due to some contraction of the flow resulting 
from the boundary layers along the walls which were 
not completely removed by the suction. Furthermore, 
this profile allowed for a comparison with theoretical 
calculations. 

A comparison of the experimental pressure distribu 
tion with theory is shown in Fig. 9 for an unstaggered 


cascade with solidity ratios dc = 1.0, 0.75, and 0.50; 
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in each case the direction of the inflow was parallel to 
the chord of the blade. The theoretical pressure dis- 
tribution over the single airfoil, d/c = 
The difference between the pressure 


©, is also shown 
for comparison. 
distribution of the blade in cascade and the single air- 
foil naturally increases with closer spacing. For all 
solidity ratios, the agreement between theory and ex- 
periment is very good indeed. 

Some of the results on loss coefficients are shown in 
Fig. 10, where the loss coefficient 


¢ = Ah/(p/2)w,” = Ah/g, (6) 


based on the velocity head of the effective flow g, = 
(p/2)w,” is plotted against the deflection Aw,/w, for 
three solidity ratios d/c = 1.0, 0.75, and 0.50. In the 
case of the unstaggered cascade, 8, = 90°, all cas- 
cades are compressor cascades. For the staggered cas- 
cades, 8, = 120°, the negative values of Aw, corre- 
spond to turbine cascades and the positive values of 


Aw, to compressor cascades. The theoretical loss co- 


AERONAUTICAL 


SCIENCES MARCH, 1954 
efficients are also shown for comparison with the ex 
perimental values. The agreement is good. Particu- 
larly, the rapid increase of the loss coefficients with in- 
creasing deflection Aw,/w, due to separation is well 
reproduced by theory. 

So far, only two-dimensional flow problems in cas- 
cades have been discussed. From the practical point 
of view, the three-dimensional losses due to the side 
walls of the passage and due to the clearance are of 
equal importance. For normal values of the ratio of 
blade height 6 to blade chord c, of b/c = 2 to 3, the 
two-dimensional losses—so-called blade losses——are 
about half the total losses. The other half is due to 
the three-dimensional effects at the walls and _ the 
clearance. Extensive results about the three-dimen- 
sional losses at the walls have been obtained earlier by 
Carter.?! 
the clearance have been given in a report by Schaffer. ' 


Some contributions about the losses due to 


Some results are given in the German version of this 


paper. 
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Fic. 11. 
number of blades, z = 44; and (d) hub, various Ry. 
radial stations r/R; blade section GHH-B. 


Experimental setup for cylindrical cascades. 


(1) General view: 
(II) View in upstream direction 


(a) inlet nozzle; (b) outer ring; (c) cascade, 
(III) Cross section of blades at various 
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Fic. 12. 
Ry/R; blade section GHH-B. R = 350 mm.; c = 63 mm.; 
constant angle of flow at outlet, Bs = 24°. 


(3) INVESTIGATIONS ON CYLINDRICAL AXIAL-FLOW 
CASCADES 
Extensive experimental investigations have been 


carried out on cylindrical axial-flow cascades in order 
to clarify the problems of three-dimensional cascade 
Scholz'* has reported on this subject in greater 
These 


flow. 
detail, and Schaffer'® has extended this work. 
investigations, extending over some years, had orig- 
inally been conducted in order to verify experimentally 
the theory of the dead flow region near the hub in a 
rotating flow as given by Bammert.'® I do not want 
to give any details about this point but want only to 
report on the relation between two-dimensional and 
three-dimensional cascade flow, which we have drawn 
from these experiments. 

The test arrangement is shown in Fig. 11. The flow 
is confined to the annular space between the hub (ra- 
dius R,,) and the outer tube (radius R). The stator 
contains 44 untwisted blades radially arranged. At 
the entry to the stator, the flow is purely axial, the 
whereas the angle of 
Thus there is a con- 
The flow leaves the 


angle of inflow being 6, = 90°, 
the outflow is about 8. = 24°. 
siderable deflection of the flow. 


Stator with a swirl of such a magnitude that, in the 
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Radial distribution of pressure difference p; — p.(r) across the cascade for cylindrical cascades of various hub atiosr 
Re = w,c/v = 6 X 104 


Theory according to E. Eckert!’ assuming 


annular space behind the stator, the circumferential 
velocity of the flow is a multiple of the axial component. 
Four different cylindrical cascades have been tested. 
The outer radius, the number of blades, and the blade 
the radius of the hub R,, and, 
follows 


profile being the same, 
hence, the hub ratio R,,/R were 


[Z): 


varied, as 
(Fig. 


R,,/R = 0.8; 0.65; 0.5; 0.4 


outer radius remained con 
stant with a value of d/c = 0.8, whereas the solidity 
ratio at the hub varied from d/c = 0.64 to 0.32. In 
the annular space behind the stator the distributions of 
of the magnitude of the velocity, 


The solidity ratio at the 


the static pressure, 
and the direction of the 
the blade spacing and the radial distance. 
sure distribution measurements over the blade surface 


measured over 
Also pres- 


velocity were 


over the whole annular space were carried out. 


In Fig. 12 the radial distribution of the pressure 
difference across the cascade is shown for the four 
different hub ratios. The pressure difference across 


the cascade, Ap = ~; — po, is considerably greater near 
the hub than near the outer wall, because the pressure 
p>» behind the cascade increases from the inner part 


toward the outer part because of the radial pressure 
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——— cylindrical cascade “+ = 0,4 
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Fic. 13. Comparison of pressure distributions for the two-dimensional and the cylindrical cascade. 


o 
Ba = 51 - B, = 90°; 


distribution, which is connected with the rotational 


flow; whereas the pressure /; in front of the cascade is 
constant over the radius. The presentation in Fig. 12 
shows the remarkable fact that the radial distribution 
of the pressure difference across the cascade is inde- 
pendent of the hub ratio. It can therefore be con- 
cluded that the rotational flow is of the same type for 
all hub ratios. For comparison, the theoretical curve 
according to Eckert and Korbacher" is included, which 
was calculated for an outflow angle of 8B. = 24° cor- 
responding to two-dimensional cascade flow. This 
theoretical curve is somewhat below the experimental 
values of the cylindrical cascade. The difference be- 
tween the two curves represents the additional pres- 
sure difference due to friction. This additional pres- 
sure difference is approximately constant along the 
radius. 

With the aid of the pressure distribution measure- 
ments it is possible to throw some light on the “radial 


divergence,’ this meaning the variation of the solidity 




















Two-dimensional cascade 


ei R 
cylindrical cascade — = 04 


cylindrical cascade 


Parameters as in Fig. 12 
Re = w,c/v = 6 X 104 


ratio along the radius. For this purpose also, pressure 
distributions on the two-dimensional cascade with the 
same blade section and with different solidity ratios 
A comparison of the pressure dis- 


and the 


have been taken. 
tributions of the 
cylindrical cascade is shown in Fig. 13. 
distributions for two values of the solidity ratio, dc = 
0.75 and 0.50, of the two-dimensional cascade are com- 
pared with those of the cylindrical cascade at such a 
radial station, for which the solidity ratio is just the 
The pres- 


cascade 
The pressure 


two-dimensional 


same as that of the two-dimensional cascade. 
sure distribution of the two-dimensional cascade is in 
excellent agreement with that of the cylindrical axial- 
flow cascade. It can therefore be concluded that the 
effect of ‘‘radial divergence’ on the pressure distribu- 
tion and, hence, on the boundary layer is negligible. 
It can further be concluded that for a stator row cas- 
cade the local loss coefficients are independent of the 
“radial divergence’ due to the absence of the centrif- 


ugal forces. 
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In Fig. 14 the local loss coefficient ¢ = Ah (p, 2) wi" 
{based on the inflow velocity) is plotted against the 
radial distance for a cascade with the hub ratio R,, R = 
0.5. Except for the region very near to the hub 
and the outer wall of the annular space, the loss co- 
efficient of the two-dimensional cascades is in full 
agreement with that of the cylindrical cascade. The 
additional end losses are confined to the vicinity of the 
inner and the outer walls. In the present case, where 
the ratio of blade height } to blade chord ¢ is b/c = 2.8, 
the total loss, including three-dimensional losses, is 
Crotat = 0.480, while the losses of the two-dimensional 
cascade are 0.391. The additional end losses there- 
fore are 19 per cent of the total losses. 

Finally, a few remarks may be made about the radial 
distribution of the velocity components behind the 
cylindrical cascade. The radial distribution of the angle 
of outflow B and the radial distribution of all three ve- 
locity components—1.e., the circumferential component 
W%., due to rotation; the axial component we,, which 
gives the mass flow; and also the radial component 
2, are shown in Fig. 15 as measured in a plane 0.03R 
behind the cascade, with a hub ratio R,,/R = 0.5. The 
radial component is scaled up by a factor of ten, and 
all quantities are averaged over one spacing in the cir- 
cumferential direction. 

For the axial and circumferential components of the 
velocity, theoretical values are given which have been 
calculated according to Eckert" for the measured dis- 
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tribution of the outlet angle 6. This theory shows 
satisfactory agreement with measurements. As _ the 
radial component of the velocity is comparatively 
small, the correct measurement of these values is rather 
difficult. 
ponent is directed toward the center throughout the 


The sign 


In the present case the radial velocity com 


entire radial distance of the annular space. 
of the radial velocity component is easy to realize 
Because of the radial pressure increase due to the ro 
tating flow (Fig. 12), there is a greater pressure differ 
ence at the inner parts of the cascade than at its outer 
parts. Consequently, the amount of the outflow ve 
locity we is greater at the inner part than at the outer 
part of the cascade. Since the direction of the outflow 
velocity ws is almost constant over the radial distance 
due to the untwisted blades, the axial component w», 
giving the mass flow is also greater at the inner part 
than at the outer part. As the axial component in 
front of the cascade is constant over the radius, a 
transportation of fluid must take place inside the cas 
cade from the outer part toward the inner part. This 
transportation of fluid is associated with an inward 
radial component of the velocity inside the cascade, 
which is still persisting shortly behind the cascade, as 
can be seen from our experiments. <A theoretical curve 
of the radial velocity distribution is also given which 
is based on the concept of the flow just discussed and 
which was calculated by applying the equation of con- 
tinuity to the w», distribution. 
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(4) CONCLUSIONS 


The results of the investigations of the flow through 
cascades are now summarized, and a few suggestions 
are given for further research in this particular field. 

Systematic investigations of cascade flow carried out 
at the Institute of Fluid Mechanics at the Technical 
University of Braunschweig were primarily done for 
two-dimensional cascades in incompressible flow  in- 
cluding viscosity effects, because we think that the 
effects of compressibility and the effects due to three- 
dimensional flow (end losses, losses due to clearance, 
and effects due to radial divergence) can be attacked 
successfully only if the two-dimensional incompressible 
cascade flow is sufficiently well understood. With re- 
gard to this problem we have succeeded in making some 
Some results have also 


remarkable contributions. 


been obtained for the three-dimensional cascade flow. 


These are the main results of our investigations: 


Two-Dimensional Cascade Flow 


(1) Theoretical.—Simple methods for solving the 
“Indirect Problem” and the ‘Direct Problem” 
been worked out for the two-dimensional incompressible 
frictionless flow through cascades. By this method the 


time for calculating the pressure distribution for a given 


have 


distance R 
! “2a 
0,7 0,8 A '~ plane of 


cascade, for ten different angles of inflow, has been 
reduced to about 20 hours. 

(2) Theoretical.—Extensive systematic investiga- 
tions of the loss coefficients have been carried out for 
unstaggered cascades in two-dimensional incompres- 
sible viscous flow. The theoretical calculation of the 
loss coefficients leads to cascades of best efficiency 
1.e., cascades with a minimum loss coefficient at a given 
angle of deflection. These cascades in most cases have 
a comparatively large value of the solidity ratio and 
are associated with a flow that is slightly separated. 

(3) Theoretical.—By applying the Prandtl-Glauert 
rule, a method has been developed for calculating the 
compressible frictionless two-dimensional flow through 
cascades with small angles of deflection of the flow 
through cascades. It has been shown by some ex- 
amples that at a constant blade section and constant 
solidity ratio the velocities at the blade in cascade are 
increasing more rapidly with increasing Mach Number 
than the velocities on the single airfoil. 

(4) Experimental.—A two-dimensional cascade tun 
nel for incompressible flow has been constructed pri- 
marily for checking the results of theoretical cascade 
calculations, but it also was constructed for investi- 
gating special cascade problems of practical impor- 


tance. Special care has been taken to obtain conditions 
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INVESTIGATIONS 
of two-dimensional flow. By means of boundary-layer 
slots and by applying suction, this has been achieved 
satisfactorily. Measurements in this cascade tunnel 
are in good agreement with theoretical results. 


Cylindrical Axial-Flow Cascades 


(5) Radial Divergence.—The radial divergence of the 
blades has a very small influence on the flow through 
evlindrical cascades. The pressure distribution over 
a blade in a cylindrical cascade was found to be exactly 
the same as in the two-dimensional cascade of the same 
solidity ratio. 

6) Radial Velocity Components. 
cascade with untwisted blades and a purely axial in- 


For a cylindrical 


flow, there exists a radial velocity component immedi- 
ately behind the cascade directed toward the center 
with a maximum value of about 25 per cent of the axial 


velocity component for hub ratio R,/ R = 0.5. 


Further Investigations 


For further investigations in the field of cascade flow 
the following problems seem to be the most important. 

(1) Experimental Investigation of Compressible Flow 
Through Cascades.—Some careful two-dimensional com- 
pressible cascade experiments are needed as the basis 
for theoretical work in this field. These measurements 
should start from Mach Numbers in the incompressible 
flow range and go up to supersonic flow, thus giving a 
link between incompressible and compressible two- 
dimensional cascade flow. 

(2) Theoretical Investigations of Compressible Flow 
Through Cascades.—It is desirable to enlarge consider- 
ably upon the calculations of the compressible fric- 
tionless two-dimensional flow through cascades, first 
in the range. The Prandtl-Glauert rule 
applies to purely subsonic flow and to moderate angles 
However, for turbine 


subsonic 


of deflection of the flow only. 
cascades large angles of deflection are required. 
(3) Three-Dimensional Losses..-Apart from the 
losses in two-dimensional cascades, a considerable con- 
tribution to the total losses results from the three- 
dimensional losses that are due to the clearance and to 
the sidewalls of the blade passage. Having analyzed 
the two-dimensional losses, a systematic investigation 
of the three-dimensional losses is urgently required. 
(4) Investigations on Rotating Cascades.—Experi- 
ments with rotating cascades are of prominent impor- 
tance, first in the region of subsonic flow. These ex- 
periments are necessary for the investigation of the 
effects of centrifugal forces on the flow through the 


cascade 
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Applications of Hypersonic Small-Disturbance 
Theory 


MILTON D. Vax DYKE* 
Ames Aeronautical Laboratory, NACA 


SUMMARY 


When slender objects fly at high supersonic speeds, linearized 
theory breaks down because nonlinearity becomes an essential 
feature of the flow. A study is made of the simplest inviscid 
theory retaining the essential nonlinearity, which is that under 
The 


iccuracy is greater than in the transonic or supersonic small 


lying the hypersonic similarity rule of Tsien and Hayes 


disturbance theories and can be further increased by interpreting 
results in conformity with Newtonian impact theory 

Because of a coincidence, this hypersonic small-disturbance 
theory can be reinterpreted so as to yield a first approximation 
also in the range of linearized theory. Thus a single theory and 
a single similarity rule cover the entire supersonic range from just 
ibove the transonic zone to infinite Mach Number 

To illustrate that the theory can serve both as a simplified 
model of the full equations and as a practical approximation, 
several problems relating to bodies of revolution are solved. An 
inomaly in an earlier solution for slender cones is found to result 
from computational errors rather than any shortcoming of the 
theory. The initial flow gradients at the nose of an ogive are de 
troublesome singularity encountered in the 
Light 


termined, and a 
full problem by other investigators is found to be illusory 
hill’s technique of rendering solutions uniformly valid is then used 


to calculate the third term in the power series for surface pressure 


SYMBOLS 


b = ratio of cone slope to shock-wave slope 


B = function defining body shape 

C, = pressure coefficient 

f = conical stream function for cone 

g,h = subsequent terms in series expansion of stream 
function for ogive 

V/ = free-stream Mach Number 

p = pressure 

= radial coordinate 

R(x = radius of meridian curve of ogive ot revolution 

s(x = function defining shock wave on body of revolution 

= = function defining shock-wave surfaces 

( « = Cartesian velocity components with wu in free 
stream direction 

be Hs = Cartesian coordinates with x in free-stream direc- 
tion 

8 = Vi? - 1 

1 = adiabatic exponent of gas (1.405 in examples) 

6 = local slope of two-dimensional body 

6 = conical coordinate, 7/* 

p = density 

o = semivertex angle of cone 


T = thickness parameter of body; in examples, initial 
slope of shock wave 


= stream function for body of revolution 


Presented at the Aerodynamics Session, Annual Summer Meet 
ing, IAS, Los Angeles, July 15-17, 1953. Revised and received 
November 26, 1953 

* Aeronautical Research Scientist 


w = reduced entropy function, p/p” 
= free-stream value 
( = value at tip of ogive 


reduced form |see Eq. (2 


ll 


jump in quantity through shock wave 


INTRODUCTION 


eres SHAPES ARE ordinarily most efficient 
when the flow disturbance. 


For this reason, linearized theory, based on the assump 


they cause least 
tion of small disturbances due to thin or slender bodies, 
has been of great practical value in analyzing both sub 
sonic and supersonic flows. However, when the flight 
Mach Number of even a very slender object greatly 
exceeds unity, linearized theory breaks down because 
nonlinearity becomes an essential feature of the flow. 
At such speeds, just as in the transonic range, linearized 
theory must be abandoned in favor of an approximation 
that retains the essential nonlinearity. 

If viscosity can be neglected (or accounted for sep 
arately), the velocity disturbances produced by a thin 
body remain small compared with its flight speed no 
matter how great the Mach Number. The pressure 
disturbances are also small compared with the stagna 
tion pressure, though perhaps large compared with the 
In this sense, the flow field involves 
Consequently, the full equa- 


static pressure. 
only small disturbances. 
tions of motion can be simplified to provide a useful 
approximation at arbitrarily high flight speeds. 

An investigation is made here of the simplest ap- 
proximation to the full equations which retains the non 
linearity essential at Mach Numbers large compared 
with unity. The resulting theory is that which under 
lies the hypersonic similarity rule of Tsien’ and Hayes. 
Although the similarity rule has already been thor 
oughly explored, the underlying theory has, in fact, 
never been written down except for plane flow in a 
restricted range of Mach Number,* and applications 
have been limited to the few special cases that are 
strictly irrotational or can be assumed approximately 
so. Here other examples are studied in which flow ro 
tation is important. It is shown that the small-dis 
turbance theory serves as a helpful guide to the full 
theory and as a useful approximation for practical cal 
culations. 

A more detailed exposition of some parts of this paper 


is given in reference 4. 
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HYPERSONIC SMALL- DISTURBANCE THEORY 


Consider a thin or slender three-dimensional body 
fixed in a uniform supersonic stream. The motion will 
be assumed steady for simplicity, the extension to 
slightly unsteady flow being straightforward.‘ Let the 
thickness of the body be described by the parameter 7, 
which is small compared with unity. (For example, 7 
may be taken as the maximum streamwise slope of the 
body.) It is known that linearized theory yields an 
adequate description of the flow only when the free- 
stream Mach Number J/ is some moderate multiple 
of unity, so that W7r < 1. 


\Jr approaches and exceeds unity, linearized theory 


As M increases, so that 


breaks down in general (and so do higher approxima- 
tions based upon it, such as second-order theory). 
Thus nonlinearity becomes essential, and the flow will 
be termed hypersonic if the free-stream Mach Number is 
so great that 1/7 is not small. (It may be very large, 
as when lf —> o~.) 

Mathematically, it is convenient to regard the vari- 
ous small-disturbance theories as approximations that 
vanishing From this 
point of view, we seek a simplified theory that provides 


become exact for thickness. 


a first approximation as 


—> ( , 
: , l with 1/Mr = 0(1)* 


M— of 


Small-Disturbance Problem 


A first approximation to the full equations of motion 
is to be found by retaining only leading terms in the 
thickness parameter 7. This is accomplished by intro- 
ducing new variables that are of order unity. 

The full equations of motion, expressing the conserva- 
tion of mass, momentum, and energy, may be written 


as 
(pu), + (pv), + (pw). = 0 (la) 

uu, + vu, + wu, + (p,/p) = 0 (1b) 

uv, + vw, + we, + (p,/p) = 0 (Ie) 

uw, + vw, + ww, + (p./p) = 0 (1d) 


u(p/p’), + u(p/p”)y + w(p/p’), = 0 (le) 


where subscripts indicate differentiation. 

An appropriate transformation of variables is sug- 
gested by the known facts that lateral velocities must 
be proportional to the body slope, that at high Mach 
Numbers the lateral extent of the disturbance field is a 
moderate multiple of the body thickness, that according 
to Newtonian impact theory the surface pressure co- 
efficient is proportional to the square of the inclination, 
etc. Thus, choosing the x axis in the free-stream direc- 
tion, we introduce new (barred) variables according to 


* Here the order symbols are used in the conventional sense: 
f(r) = O(1) as r > O means that f(r) remains bounded as rt ~ 0, 
and f(r) = O(g[7]) means that f(7r)/g(r) is O(1). 
f(7) = o(1) means that f(7) vanishes as r > 0. 


Similarly, 
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u=u,[1 + 7r°a(%, J, Z)] 
v = u,td(2, Y, 2) 
w= uTws, Y, Z) 


b = p,(yM?r’)p(, 9, 2) 
Pp = p.Ap(&, ¥, Z) 
r 


We assume that the new dependent variables (i, 
i, @, p, p) are O(1) as r > O for fixed Mr and that the 
reciprocals of the new independent variables (1 J, | 7, 
1/Z) are alsoO(1). This assumption will be justified by 
the complete internal consistency of the resulting ap- 
proximate problem. 

If we introduce these new variables and discard terms 
that contain 7° explicitly, Eqs. (1) reduce to 


a; + da, + @iz + (p/p) = O 3 
and 
Dz + (p7); + (PB); = O (4a 
d; + i, + @i; + (p;/p) = 0 tb 
@, + o@, + @®; + (p/p) = 0 (4c 
(P/B)s + H(b/B’)s + @(b/B")3'= O (4d 


Upstream and laterally far away from the body, the 
flow quantities must attain free-stream values, 


uw=v=%H=0 (ja) 
b = 1/yM?r in the free stream (5b) 
p = 0 (5c) 


If the body is described in reduced coordinates by 
B(z, ¥, 2) = 0, the condition of tangent flow at the sur- 
face becomes, neglecting terms of 0(7”), 

B,.+%t7B,+@#B,;=0 (atB=0 (6) 

If the full system of shock waves is described by 
S(z, 3, 3) = 0, the shock-wave relations become, to the 
same approximation, 

a] Ss = |7 bY = [w] S; (7a) 


[p(.S; + aS; + WS;) | = 0 (7b 


(a(S; + dS; + @S;)? + (S,? + S2)p] = 0 (Te) 


i; (S; + aS; + @S;)? + (S;° + S;*) p | = () 
2 7 

(7d) 

[p/p] > 0 (7e 


at S = 0. 
in a quantity across the shock wave. 


Here square brackets denote the increase 


The assumptions regarding orders of magnitude are 
now justified by the complete internal consistency of the 
reduced problem when | (1/7) = O(1). However, the 
consistency argument breaks down when J/7 < | 


([formally, when \f7 = o0(1)], because then the require- 
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1 in the free stream 
This is an 


ment of Eq. (5b) that p = (7.1/*r") 
contradicts the assumption that p = O(1). 
automatic warning that the reduction breaks down in 
the range of linearized theory. 

A significant feature of the reduced equations is that 
the cross-flow problem for (@, @, p, p) does not involve 
the streamwise velocity #7. Hence the cross-flow prob- 
lem can be solved independently; then wu (if required) 
is given by the reduced Bernoulli equation, 


l ; a Y Pp 
6 Pa i ) ~ = const. (S) 
- 7 l p 
Therefore Eq. (3) and the first terms of Eqs. (5a) 


and (7a) are superfluous and can be disregarded hence 
forth. 

The cross-flow problem is completely equivalent to 
the full problem of unsteady motion in the j-2 plane 
due to a moving solid boundary described by (J, 
This analogy was pointed 


j, 2) = O, # being the time.* 
out by Hayes in his extension of the Itypersonic simi 
larity rule.* 

Because terms of order r* have been disregarded in 
the reduction, flow quantities determined from the 
small-disturbance equations will differ from their full 
values by O(7*). This more general result agrees 
with Goldsworthy’s conclusion® that in the special case 
of plane flow with 1/7 ~ 1 the error is O(1/.1/*). The 
present estimate applies, however, even at lJ = o. 


Combined Supersonic-Hypersonic Small-Disturbance 
Theory 
It is known that transonic small-disturbance theory 
embraces linearized subsonic and supersonic theory as 
a special case, so that the awkward question of de- 
limiting the transonic zone need not arise.> Likewise, 
a combination of hypersonic small-disturbance theory 


* The undisturbed fluid having density px = 1 and pressure 


Do = (yM?r?)—! 
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with linearized supersonic theory can be effected® but 
only through a coincidence. An immediate obstacle to 
such combination is the fact that the approximations 
leading to the two theories are mutually exclusive 
The difficulty arises in the continuity equation, Eq. 
(la), which in hypersonic theory reduces to 


UaoPr + (pv), + (pw), = O 
but in linearized theory reduces instead to 


UoPr t+ p.(Ur + Vy + we) = O 


The term p..u,; must be retained in linearized theory, 
whereas it must be dropped in hypersonic theory in 
order to achieve similitude. 

However, comparison of the similarity rules in the 
two cases indicates that the two theories may be arti 
ficially combined simply by reinterpreting the solutions 
According to the hyper 
sonic similarity rule, the reduced variables of Eq. (2 
depend upon .\/ and 7 only in the combination 1/7. 
It can be shown that solutions of hypersonic small- 
disturbance theory remain valid for small values of 
\{r (which is the domain of linearized theory), pro 
vided that \/7r is replaced by 87 and the results are 
interpreted in terms of physical variables according to 


of the hypersonic theory. 


u = u.(1 + 7r7u(%, J, 2; Br)] 
v = UuTi(«, ¥, 2; Br) 
w= uU,Tw(r, ¥, 2; Br) 
pb = pol 1 + yM*?? BLE, 9, 2; it) = 

1 ) (9) 

- =p M?r°p — 
aT" | ‘ Y t : 

q M? , — , 

p = Ppa Ll + 3 D(z, J, 3; Br) —13 |= 





rather than according to Eq. (2). This means, for ex- 
ample, that the hypersonic similarity rule for C,, 7° 
covers the entire range of supersonic flows if Wr is re 


placed by @r. 
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The error in hypersonic small-disturbance theory is 
O(7r*), and that in linearized theory is 0(7/8) for shapes 
such as flat wings but only 0(7°/8*) for slender objects 
such as bodies of revolution. Hence the error in the 
combined theory is O(7°, 7/8)* for wings and O(7°, 
7°/B*) for fusiform shapes. 

The obvious advantage of the combined theory and 
its associated similarity rule is demonstrated in Fig. 1. 
The pressure on cones of various thicknesses’ is corre- 
lated first according to the original hypersonic similar- 
ity rule and then according to the combined rule, which 
is seen to extend the range of correlation down to the 


transonic zone. 


Optimum Use of Hypersonic Small-Disturbance Theory 


The error in hypersonic small-disturbance theory is 
only O(r*) compared with O(7) for linearized theory and 
0(7r’*) for transonic small-disturbance theory. Hence 
the need to correct for higher order terms in thickness 
is less in the hypersonic range than at lower speeds. 
Nevertheless, the accuracy of hypersonic small-dis- 
turbance theory can be strikingly improved simply by 
interpreting its results in conformity with Newtonian 
impact theory. 

Newtonian impact theory assumes that fluid par- 
ticles striking the body lose the normal component of 
their free-stream momentum, with the result that sur- 
face pressures are proportional to the square of the sine 
of the streamwise inclination. Although this result 
is not exact even in the limit as JJ ~ o (and y > 1), 
it holds approximately at high supersonic speeds. 
This suggests that the pressures predicted by hyper- 
sonic small-disturbance theory be extrapolated to 
thicker shapes by assuming a sine-squared variation. 

Fig. 2 shows that in the limit of 1J = o this idea ex- 
tends the utility of the hypersonic correlation to cones 
of semivertex angles as extreme as 50°. Another illus- 
tration arises in correlating the lift on inclined bodies 
of revolution. Hypersonic small-disturbance theory 


* Whichever is the greater. 














Cp 
sinto 
1.6 
Cp 
tanto 
8 FT T 
O 
O 10 20 30 40 50 
o, deg. 
Fic. 2. Correlation of pressure on thick cones at MM = o@, 


THE AERONAUTICAL SCIENCES 


MARCH, 1954 


2.2 


Normal 
force 


“Wy 


Modified 
component 








2.0 
® 
a 
| 
” 
@ 
4 
R= 
— hey 
o 
F 
S 
2 
1.6 
o 667 | 2 @ 
B tane 
Fic. 38a. Correlation of lift of cones. Normal force 
€ 2.2, 
r= 
6 ° 
a 
Ee 
3 
rs) 
a] 
a4 
= 
° 
E 
a) 
© 
a 
2 
” 
1.8 : 
x .667 | 2 @ 
B tana 
Fic. 3b. Correlation of lift of cones. Modified component 


cannot distinguish between lift and normal force. 
However, assuming a sine-squared variation of pressure 
indicates that, whereas the initial lift and normal- 
force slopes both vary considerably with body thick- 
ness, the modified component that is rotated twice as 
far behind the lift vector as is the normal-force vector 
(see sketch on Fig. 3a) should be independent of thick- 
ness. Fig. 3 shows that for cones’ this is so nearly true 
that the accuracy of correlating hypersonic lift can be 
increased many fold. 


APPLICATIONS TO BODIES OF REVOLUTION 


The utility of hypersonic small-disturbance theory 
will now be illustrated by applying it to three related 
problems connected with pointed bodies of revolution. 
The first of these problems has been completely solved 
for the full equations, so that comparison will indicate 
the accuracy to be expected from small-disturbance 
theory; in the second problem the small-disturbance 
theory has recently served as a guide to the full solu- 
tion; and in the third problem useful results are obtained 
which have not yet been attempted with the full theory. 

For axially symmetric flows it is convenient to satisfy 
the reduced continuity equation [Eq. (4a)] by intro- 
ducing a stream function y, 
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vi = 7p (10 
y= 


{ 
pol 
Then Eq. (4d) states that p p’ is a function only of 
y (that is, entropy is constant along stream lines be- 
tween shock waves). Substituting into the momentum 


equations gives 


Yi Wes — 2WibiWer + Via Wri 
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where 
w(y) = p p’ (11b) 


The problems to be considered involve a single bow 
shock wave, which may be described in reduced co- 
ordinates by ? = s(#). Then the shock-wave relations 
obtained from Eqs. (7) (or simply by reducing the con- 
venient results of reference 9) are 

2 M?*r’s’? — 1 
7 = (12a) 


¥+ 1 M?r?5’ 


2yMer?s’? — (y — 1) : 
p= -atr = s(t) (12b) 


y(y¥ + 1) AP?? 


(y + 1)AM?r?s’? 
p= (12c) 


2+ (y — 1).M?*r’s” 


Circular Cone 


Consider flow past a slender circular cone of semi- 
vertex angle o. It is convenient to regard the Mach 
Number and shock-wave angle as given and to identify 
the thickness parameter 7 with the slope of the shock 
wave, letting ¢o = br, where the constant 6 must be 
found. 

The flow field is conical, so that the stream function 
has the form 


¥(Z, 7) = <°f(0) (13a) 
where 
6 = </? (13b) 


is a conical variable varying from 6 on the cone surface 
to unity at the shock wave. Then Eq. (lla) becomes 
the nonlinear ordinary differential equation, 


we aA LAU f’ ofits . cw f' 
4f?f — Off"? = yoo = if" = - (14a) 
where wy is the constant value of £/p” behind the shock 
wave given, according to Eqs. (12b) and (12c), by 


2yM?r? — (y — 1) f2 + (y — IM??? Pp 
ae es (14b) 
y(y + 1) M?*?? (y + 1) M?r? 


The boundary conditions at the shock wave are 
found, from Eqs. (12a) and (12c), to be 
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(y + 1)M?r? l4e 
2+ (7 -— 1) ee} 


The constant 6, which gives the cone surface, is deter 
mined by the requirement that the stream function 


vanish on the surface, 
f(b == {) (l4d 


The nonlinear equation, Eq. (14a), is readily inte 
grated numerically, starting at the shock wave (@ = | 
and proceeding inward until f vanishes. The pressure 


coefficient is given in terms of f’ by 


C wy A) | : 
= 9 aa (lo 
a? | 4 y Mo" 


In this way the surface pressure has been determined for 
values of the similarity parameter greater than 
(Below this range, linearized theory becomes accurate. 
These results, expressed in the combined supersonic 
hypersonic form, are compared in Fig. 4 with the full 
solution’ for cones of various thicknesses. The frac 
tional differences between the approximate and full 
solutions are nearly equal to 7°, which suggests that 
the estimate of the error as 0(7*, 7*/8") can be relied 
upon quantitatively. (The differences could, of course, 
be considerably reduced by interpreting the results in 
conformity with Newtonian impact theory. ) 

This approximate solution for cones was previously 
sought by Shen,'® whose result is also shown in Fig. 4. 
It appears that his solution must involve computational 


errors. 
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Initial Gradients on Ogives 


The initial gradients of flow quantities at the tip of 
an ogival body of revolution can be determined by per- 
turbing the solution for a cone. It is clear that the 
stream function will then consist of that for the cone 
plus an additional term, 


W(x, 7) = r°7f(0) + 72(8) (16) 


Substituting into the equation of motion, Eq. (lla), 
yields a linear ordinary differential equation for the 
perturbation function g, with nonconstant coefficients 
depending upon the basic function f. Specifying the 
initial curvature of the shock wave provides initial 
conditions upon g. The equation can then be inte- 
grated numerically, and thus the initial pressure gra- 
dient and the ratio of shock wave to body curvature 
are determined. 

A complication arises near the body surface, because 
6 = bisa singular point of the differential equation for 
g. Asa consequence, it is found that g is not analytic 
near the surface but includes a term proportional to 
(@ — bd) This term arises from the fact that the 
pencil of fluid striking the tip of the ogive is spread thin 
over the entire surface, and the linear entropy gradient 
at the tip due to a curved shock wave is thus intensified 
to a square-root gradient normal to the surface else- 
where. 

This problem has recently been considered for the 
full equations by Cabannes!' and by Shen and Lin." 
The exceptional behavior of the solution near the sur- 
face was glossed over by Cabannes. On the other 
hand, Shen and Lin claimed to have found a logarith- 
mic singularity. In addition to other difficulties, this 
would imply that the initial pressure gradient is in- 
finite, although numerical solutions by the method of 
characteristics give no such indication. If the present 
small-disturbance theory is a faithful model of the full 
equations, the absence from it of any logarithmic singu- 
larity means that none can actually appear in the full 


| ] Fai 
N Cone-expansion approx. / 

: J 
hs 

te id 


Small-disturbance theory 


6 N 














4. } 
aC, 

a(R,Rex) r re Rox + eRox® + ERO’ x? +... 
e Ff x 
(@) lies 

RK Fo) a | ie@) 
BRo 
Fic. 5. Initial surface pressure gradients on ogives. 


THE AERONAUTICAL SCIENCES 


MARCH, 1954 








15 
| 
| 

| 

| 

















10 | 
a*Cp 
we 
(RS x) | 
Small-disturbance theory 
5 \ 
| 
ye | 
Cone-expansion approx. 
| | 
| 
Peo Fe o @ 
BR, 
Fic. 6a. Initial surface pressure curvatures on ogives 
N Cone-expansion approx. / 
i 
Small-disturbance theory 
4 4 
ap 


B(R, x) (Rex) 





2 } 
ie) 
.333 co) I © 
‘ 
BR, 
Fic. 6b. Initial surface pressure curvatures on ogives 
problem. This conclusion has been confirmed by the 


authors of reference 12, who find upon re-examination 
that the logarithmic singularity was in fact illusory.* 

The initial pressure gradient is shown in Fig. 5 for 
values of the combined similarity parameter (based 
on initial slope) greater than ! 


Further Terms in Series for Ogives 

It might be supposed that, for an ogive described by 
a power series, further terms in a power series expansion 
for surface pressure could be found by retaining addi- 
tional terms in Eq. (16), 


¥(z, 7?) = £°f(0) + z*g(0) + F4h(0) +... (17 


However, complications arise at this stage because of 


* Private communication; see also Addendum No. 1 to refer 


ence 12 
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the nonanalytic behavior of the stream function near 
the surface. Whereas g(@) contains an authentic */>- 
power branch point, /(@) is found to contain a spurious 
power branch point, the next term an inverse !/>- 
power singularity, and so on. Hence, straightforward 
continuation of the procedure breaks down. 

The difficulty arises from the fact that in the first 
perturbation the branch point in g arises at the basic 
cone (6 = 6) rather than at its proper location on the 
ogive surface. Although this discrepancy is _ negli- 
gible in the first perturbation, it is compounded in sub- 
sequent terms so as to be catastrophic. The remedy 
is to choose a slightly curvilinear coordinate system 
such that for each term the */s-power branch point 
Both the diffi- 
culty and the remedy are just those considered by Light- 
hill in his discussion of a technique for rendering ap- 


appears precisely at the body surface. 


proximate solutions uniformly valid.'* ' 

With this modification, the procedure can, in prin- 
ciple, be continued indefinitely. It can therefore be 
concluded that an analytic body of revolution is ac- 
companied by an analytic bow shock wave at super- 
sonic speeds and an analytic pressure distribution. 

Calculation of the next term in the series for surface 
pressure has been carried out,‘ and the results are pre- 
sented in Fig. 6. Two functions arise, the surface 
pressure being given by 


C, = Ry’? Cy + Ry’Ry”’ oe 
— oe La 0 La(Ro’Ro'’x) J. 
US pyre] OS | pip x 
2 Lorine) 
rt, 


+... (18) 


O(Ro’x)O(Ro'’’x) 4 


where Ry’, Ro’’, and Ro’’’ are the actual slope, second 
derivative, and third derivative at the tip of the ogive. 

Newtonian impact theory, which was seen to be a 
useful qualitative guide, suggests that this number of 
terms suffices to predict the entire pressure distribu- 
tion over circular-are ogives. Accordingly, Rossow’s 
accurate numerical solutions for circular-are ogives™ 
at a similarity parameter of unity are compared in 
Fig. 7 with the small-disturbance series. Three terms 
are seen to describe the pressure distribution with good 
accuracy over less than half the ogive length; at least 
one more term would be required to give the entire 


distribution. 


CONCLUDING REMARKS 


It has been seen that hypersonic small-disturbance 
theory serves both as a simple guide to the full theory 
and as a useful approximation for obtaining practical 
results. It appears that any extensive solution of 
supersonic flow using the full equations should be pre- 
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ceded by a pilot solution according to small-disturbance 
theory. 

The small-disturbance equations are still inherently 
difficult to solve, because they are essentially non- 
linear. Further approximations are desirable for sim 
plifying the theory. 
tremely simple and accurate approximations are avail 


In the case of plane flow, two ex 


able: the reduced shock-expansion theory,'® ' and 
the even simpler tangent-wedge approximation which, 
when written in combined supersonic-hypersonic form, 


gives 


! 


+1 /y +1\? 
Cc, = 6° ; + VC ) — ; (19) 
2 2 3°6° 


if 6 is the local slope (cf. reference 18). Hence, in plane 
flow the inviscid hypersonic problem is solved for prac 
tical purposes. 

For general three-dimensional shapes, a promising 
approximation of comparable simplicity has been ad 
vanced by Eggers and his colleagues’ '* *° and inde 
pendently (for bodies of revolution) by Zienkiewicz and 
Bolton-Shaw.*! ** This ‘“‘cone-expansion”’ method con- 
sists in applying the Prandtl-Meyer relations behind 
the tip of a pointed body as if the flow were plane. 
Thus, for a body of revolution, the surface pressure 
and Mach Number at the tip are taken from the solution 
for a cone, and the subsequent pressure distribution is 
calculated according to a Prandtl-Meyer expansion. 
This approximation yields surprising accuracy at 
sufficiently high Mach Numbers, and Eggers has dis- 
cussed the reasons for this.*? Zienkiewicz and Bolton- 
Shaw extend the approximation to lower Mach Numbers 
by using an empirically determined fraction of the full 
planar expansion. As an indication of the accuracy 
obtainable, the dashed curves in Figs. 6 and 5 show 
the coefficients calculated according to the cone-ex- 
pansion procedure (with the full planar expansion). 
The full capabilities and limitations of this approxi 
mation have yet to be determined. 
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Conference on Supersonic Flow 


An unclassified conference on Supersonic Flow will be held in the Mathematics Build- 
ing, University of Maryland, from 9:30 a.m. to 5:00 p.m. on March 19-20, 1954. 

Information concerning the conference may be obtained from Dr. E. L. Resler, Associ- 
ate Research Professor, Institute of Fluid Dynamics and Applied Mathematics, University 
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Knergy Approach to the General Aircraft 
Performance Problem 


EDWARD S. RUTOWSKI* 
Douglas Aircraft Company, Inc. 


ABSTRAC1 


The general aircraft performance problem is considered from 
the point of view of the balance that must exist between the po 
tential and kinetic energy change of the aircraft, the energy dis 
sipated against the drag, and the energy derived from the fuel 

This approach yields as one result a basic equation for the rate 
of change of the sum of the potential and kinetic energies of an 
iircraft. The form of the equation points to the use of aircraft 
total energy rather than altitude as the significant independent 
variable in the climb performance of high-speed aircraft. Using 
this equation, a method is outlined which permits finding either 
the path of minimum time or minimum fuel to change from one 
combination of speed and altitude to another 

The energy balance approach yields as another result a form 
of the aircraft range equation which incorporates a correction for 
the total energy change during the cruising portion of the flight 
This form of the range equation shows that the dimension and 
the order of magnitude of aircraft range are given by the heat 
content of the fuel. Finally, the parameters of the range equa 
tion are examined to obtain an indication of the range capabil 
ities of aircraft at supersonic speeds relative to those at subsonic 


speeds 


INTRODUCTION 


enema NEW PROBLEMS HAVE ARISEN in the per- 
formance analysis of high-speed aircraft which 
cannot easily be solved by textbook or the conventional 
performance methods. One of the more important 
of these problems is the determination of the optimum 
path for an aircraft to reach speed and altitude, where 
optimum may imply either minimum time or minimum 
fuel. For low-speed aircraft, it was possible to state 
the problem simply as one of achieving a given alti- 
tude. Now, however, particularly with the wide range 
of speed of supersonic aircraft, the speed change in 
getting to altitude is significant relative to the altitude 
change so that the problem is how to attain both speed 
and altitude. The conventional methods of aircraft 
climb analysis using altitude as the independent vari- 
able are adequate for the problem of simply getting 
to an altitude but are not well suited for the more gen- 
eral problem. A method using an independent vari- 
able that combined speed and altitude would be better 
suited, since the problem is how to attain a combination 
of the two. Since speed can just as well be expressed 
by a level of kinetic energy and altitude by a level of 

Presented at the Aerodynamics Session, Annual Summer Meet- 
ing, IAS, Los Angeles, July 15-17, 1953. Revised and received 
September 14, 1953 

* Aerodynamics Engineer, Santa Monica Division. 


potential energy, one means of describing a combi 
nation of speed and altitude is suggested. Represented 
as energies, speed and altitude can then essentially be 
added so that the problem can be stated as one of at 
taining a level of the sum of a kinetic and a potential 
energy, which can be defined as the total energy of the 
aircraft. The total energy would therefore seem to be 
a significant quantity in the climb analysis of an air 
craft and suggested an approach to the problem from 
an energy point of view. This approach was tried and 
found to give a relatively simple solution to the prob 
lem of how to attain speed and altitude. Furthermore, 
it was found that, not only was this approach fruitful 
in problems of climb, but it also had application to prob 
lems of aircraft range as the result of its better exhibiting 
the basic elements that determine aircraft performance. 

The purpose of this paper is to describe this energy 
approach, at first in general and, then, its application to 


problems of climb and range. 


SYMBOLS 


Cy, Ci, C2, Cs = arbitrary constants 

Cc) = airplane lift coefficient 

D = airplane drag, lbs 

E = airplane total energy, ft.lbs 

E/U = airplane specific energy, ft.lbs. per Ib 
fod te ,/, = arbitrary functions 


= acceleration due to gravity, ft. per sec 


h = altitude, ft 

H = heat content of the fuel, B.t.u. per lb 
iL = airplane lift, lbs 

V = Mach Number 

p = atmospheric pressure, lbs. per sq.ft 
R = airplane range, ft 

R = gas constant for air, ft.2/sec.2°R 

S = airplane wing area, sq.ft 

S.F.C = specific fuel consumption, lbs. /sec.lb 
t = time, sec 

7 = engine net thrust, lbs 

7 = atmospheric temperature, °R 

J = airplane velocity, ft. per sec 

J = indicated air speed, ft. per sec 

We = airplane gross weight, lbs 

W = initial airplane gross weight, Ibs 

W’ = final airplane gross weight, lbs 

W = fuel weight, lbs 

a = arbitrary variables 


: = ratio of specific heats 
n = overall powerplant efficiency 


p = atmospheric density, slugs per cu.ft. 
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THE ENERGY BALANCE RELATIONSHIP 


The energy approach to the general performance 
problem logically starts with a consideration of the 
source of the energy for propulsion. This source is the 
power-plant fuel, which is characterized by a heat con- 
tent, //7,, which is essentially an energy content per 
pound of fuel. The stored energy of the fuel is con- 
verted into work of propulsion by the power plant at 
some overall efficiency, mo, which is the product of its 
thermal and propulsive efficiencics. The energy per 
pound of fuel delivered by the power plant—that is, 
the product of the overall power-plant efficiency and 
the heat content——is expended in dissipation against 
the drag and in increasing the total energy of the air- 
craft. The sum of its kinetic and potential energies 
makes up the total energy of the aircraft and may be 


written as 


: WV? : Vy? E : 
k= Wh+ = Wt{h+ = -IW (1) 
2¢g 2¢ i 


where 
(E/W) = h + (V?/2g) 


The application of the principle of conservation of 
energy over a time interval dt, yields the following equa- 


tion: 
dE = nlldW, — (E/W)dW, — DV dt 


which describes the basic energy balance that must 
exist between the total energy change of the aircraft, 
the energy derived from the fuel, the total energy 
change due to a weight change in the airplane, and the 
energy dissipated against the drag. The terms in- 
volving the change in the total energy of the aircraft 
can be combined by use of the equation obtained by 
the differentiation of Eq. (1). 


dE = Wd(E/W) + (£/W)dW 
or 
dE = Wd(E/W) — (E£/W)dW, 
since 
dW = —dW, 
The energy balance relationship then reduces to 
Wd(E/W) = nll dW, — DV dt (2) 


This equation is here written in terms of the change in 
energy per pound of aircraft weight or what may be 
called the specific energy of the aircraft. The specific 
energy is a function only of the speed and altitude of 
the aircraft, and it, rather than the energy itself, which 
is a function of gross weight also, is the item of pri- 
mary interest in the usual performance problems. 
It should be made clear, however, that Eq. (2) does 
properly take into account the energy change due to the 
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change in gross weight as fuel is burned. Eq. (2) is 
the basic relationship from which the remaining equa- 
tions in the paper are obtained. In fact, this general 
relationship will yield most of the equations used in the 
calculation of specific performance problems. In the 
present paper, it is used to obtain the classical form of 
the climb and range equations, as well as alternate forms 
that are considered to be better adapted to the treat 
ment of some performance problems of high-speed air 
craft. 

This relationship is particularly useful in that it per 
mits a common approach to several performance prob- 
lems—for example, climb and range—which are usually 
treated as completely independent problems. 


THE CLASSICAL CLIMB EQUATION 


The climb equation in its classical form will now be 
obtained from the energy balance relationship by first 
deriving the more fundamental equation of the rate of 
change of specific energy of an aircraft. To accomplish 
this derivation, Eq. (2) is rewritten to obtain 


d(E/W) mlH.dW, DV 


dt W dt Ww 


Remembering that » J/7, represents the useful work of 
propulsion done per pound of fuel, this quantity may be 
expressed in terms of the power-plant thrust, 7, and 
the fuel rate as 
work work, time TV 
noll. = = = i : = = (+4) 
Ib. of fuel lb. of fuel/ time (dW, dt) 

Using this expression in Eq. (3), the following relation 
ship is obtained : 


d(E/W)/dt = (V/W) (T — D) (5) 


This equation states simply that the rate of change of 
specific energy is given by the excess power per pound 
of aircraft weight. It is to be distinguished from the 
commonly used relation, 


dE/dt = (T — D)V 


which is only approximate, since it neglects the energy 
change due to the changing weight of the aircraft. The 
full expression, using Eqs. (4) and (5), 


dE ds | (E/W) | 
=(T-D)Vj1- _ 

dt HT.{1 — (D/T) |no 

shows that this effect becomes increasingly important 
as the specific energy capability, or essentially the 
speed and altitude capability of aircraft increases rela- 
tive to the specific energy content of the fuel. The 
fact that the specific energy concept yields as simple 
an expression as Eq. (5) leaves little, if any, reason for 
neglecting the effect of changing gross weight. This 
relationship can be used to obtain the so-called ‘‘funda- 
mental’ performance equation of classical aerodynam 
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ENERGY APPROACH TO 
ics. The derivation is accomplished by first differ- 
entiating Eq. (1) to obtain 
ad(E/W)~ dh 4 dl 
dt dt g dt 


Using this expression in Eq. (5) and solving for the rate 


E ee | dV 
= (7T — D)i — (6) 
da sW g dt 


The portion of the equation in the brackets is the so- 
called ‘‘fundamental”’ performance equation. As shown 
by Eq. (6), it is based on the assumption that the climb 
In this particular 


of climb, 


speed be constant (dV/dt = 0). 
case, the ‘fundamental’ equation and Eq. (5) are 
identical, so that the excess power per pound of weight 
can be interpreted as giving either the rate of change of 
specific energy or the rate of climb at constant speed. 
Regarding it as a rate of change of specific energy, the 
excess power per pound of weight can also be interpreted 
as giving the longitudinal acceleration at constant alti- 
tude. In the early days of aviation when speeds were 
low, this approximate equation did not introduce ap- 
preciable error even in the case of a changing speed 
during a steady climb. The form of the correction can 
be seen by a rewriting of Eq. (6) to obtain 


dh yo. | l | = 
= (T — D) > : (4) 
dt W 1+ (Il g) (dV dh) 


The magnitude of this correction factor is shown in 
Fig. | for the simple cases of climb at constant indi- 
cated speed and constant Mach Number. Although 
small for low climbing speed, the correction factor indi- 
cated in the figure rapidly becomes appreciable as the 
climb speed increases. In fact, for climb at a constant 
Mach Number of 2.0, the correction factor exceeds 
2.0. The correction factor is well known, so that climb 
calculations of high climbing speed aircraft are gener- 
ally made using the complete expression given by Eq. 
(7 

Although Eq. (7) is adequate for the calculation of 
climb performance once the climb technique is known, 
its form is unwieldy for the more general problem of 
finding the optimum technique. The difficulty arises 
from the fact that the rate of climb is given by a func- 


tion of the type 
dh/dt = f(h, V, dV dh) 


The dependence of the rate of climb on d1’/dh makes it 
necessary to consider the climb as a whole in a problem 
such as minimizing the time to a given altitude. The 
common procedure of considering the rate of climb 
independently at each altitude is not strictly valid. 
The problem of minimizing time to climb with a func- 
tion of the above type is a typical calculus of vari- 
ations problem. However, the lack of a simple analyti- 
cal expression for thrust makes such an approach diffi- 
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cult. Simplification of the problem is generally ac 
complished by assuming that (dI’/dh) can be neg 
lected. A first approximation can be thus obtained 
which can be adjusted to give a path closer to the opti 
mum. However, this adjustment requires judgment 
so that some other means of simplification leading to a 
direct answer is desirable. It is here argued that Eq 
(5), using the concept of specific energy, provides the 
desired simplification with no loss in accuracy com 


pared to the use of Eq. (7). 


THE USE OF THE SPECIFIC ENERGY CONCEPT IN 
CLIMB PROBLEMS 


The rate of change of specific energy of an aircraft, 
at a fixed gross weight, is a function only of altitude and 
velocity. 


d(k/W)/dt = f(h, V) 


The use of Eq. (5) therefore provides the desired sim 
plification by the elimination of the derivative, d1", dh. 
This equation does not specify the form in which the 
energy is added but rather gives an expression for the 
rate of increase of the sum of the kinetic and potential 
energies of the aircraft per pound of weight. Unlike 
the ‘fundamental’ performance equation and Eq. (7), 
the climb technique, or essentially the proportions in 
which the kinetic and potential energies are added, 
does not have to be specified in any calculation using 
the concept of specific energy. The simplicity of Eq. 
(5) would therefore recommend its use in preference 
to Eq. (7) if specific energy could be shown to be a useful 
and a usable concept. 

It is to be noted that Eq. (5) is exactly the ‘‘funda 
mental’ performance equation except that the altitude 
is replaced by the specific energy. Because this spe- 
cific energy has the semblance of an altitude, the sug 
gestion has been made!’ that it be called an ‘“‘energy 
height.’ It has, in fact, the dimension of a length and 
is the sum of the geometric altitude and the altitude 
equivalent of the kinetic energy were it all converted 
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into potential energy. The significance of the use of 
the specific energy rather than altitude as the inde- 
pendent variable may be seen in Fig. 2. The contours 
of constant specific energy are the locus of all points 
having a given sum of kinetic and potential energies 
per pound of weight. Also shown in Fig. 2 superim- 
posed on the specific energy contours are the speed- 
present-day 
From 


altitude envelopes representative of a 
fighter and a hypothetical supersonic fighter. 
these envelopes, it is evident that the final altitude in 
itself is no longer adequate in the specification of climb 
problems, since this altitude may be achieved over a 
wide range of speeds. The problem is more often to 
achieve some combination of both speed and altitude 
rather than simply an altitude. Therefore, since any 
combination of speed and altitude determines a level of 
specific energy, the problem can be stated as one of 
achieving a final energy level from some initial value of 
specific energy. While it is true that a specific energy 
does not determine a unique combination of speed and 
altitude, the exact proportion of kinetic and potential 


energy may generally not be important, except per- 
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haps at the final condition desired. Furthermore, if it 
is important, the desired proportion can generally be 
achieved readily, since only a redistribution rather than 
an increase in energy is required. Because of this pos 
sibility of ready interchanges of kinetic and potential 
energy by means of a dive or a zoom, the concept of 
specific energy has tactical significance. For example, 
it would seem more important for an interceptor to 
have a specific energy capability greater than that of 
its bombing adversary than to have simply an altitude 
or speed advantage. Even in the case of other than 
tactical aircraft (research vehicles, for example), the 
concept of specific energy is useful in that it permits a 
rapid estimation of the maximum speed in a dive or the 
maximum altitude in a zoom. A first estimate can be 
made simply by finding the maximum specific energy 
capability of the aircraft in level flight. The other 
combinations of speed and altitude corresponding to 
this energy level then represent capabilities of the air- 
craft not attainable in level flight. On the basis of 
these arguments, it appears that specific energy is a 
useful concept. It remains only to show that it is a 
usable concept. 

The usability of the specific energy concept can per- 
haps best be shown by its application to a specific prob- 
lem. For this purpose, the problem chosen is that of 
finding the path of minimum time to change from one 
combination of speed and altitude to another. Using 
specific energy as the independent variable, the time 
can be written as 


-(E/U | E 
7 , y d ) (SN) 
J, w), (Ek W)/dt & 


T — D)V/W =f(V,h) 


Eut 


d(k/W)/dt = 


or since 


ee J a E i] E r] 
dt i! el led 


Eq. (S) then becomes 
7 (E/W)s E E 
§ = / fi ( af ») d ( ) or 
E/W) i HW 
“(E/H E E 
t= / fs ( a v) a( :) (9) 
E/W) HW W/ 


The problem can now be stated mathematically as one 
of minimizing the above integral. The simple condi- 
tion, easily verified by the calculus of variations,* that 


this integral be a minimum is 


* Due to E. W. Graham, Consultant, Douglas Aircraft Com 
pany, Inc. 
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ENERGY APPROACH TO 


ra) dt | 
—— - = 0 or 
ht Lad k HW’) E/W =const 


O | dt | 
= 0 (10) 
Ol’ dk IV’) EW =const 


or since, in general, for f(x, y) # O 
(0 Oy) [f(x, y)] = 0 


then 


(x, y)] = 0 


(0 Ov) [f 


and we may write 


P) [“ 7] 
= Q or 
Oh dt E,W =const 


re) k k ad 
ol’ dt EWU const 
or, finally, 


ra) it — 2) 
= Q or 
oh a IW EW const 


O Eg — DP) | 
3 : = 0 (11) 
Ol i Eu const 


The path of minimum time is therefore the one that 
This path can- 


ll 


satisfies this condition at every point. 
not be obtained analytically in general, but a simple 
graphical method analogous to the common method 
used with the ‘‘fundamental”’ performance equation can 
be used. The usual method can be mathematically 
described as one in which the path satisfies the condi- 


tion 
(0/OV) [((T — D)V/Wh), =const, = 0 (12) 


Graphically interpreted, Eq. (12) defines the procedure 
of selecting at a given altitude the speed for maximum 
climb from a plot of excess power per pound of weight 
versus speed. The interpretation of Eq. (11) is com- 
pletely analogous, except that in this case the excess 
power per pound of weight is plotted at a constant 
value of specific energy. Here, either speed or altitude 
can be used as the independent variable for purposes of 
plotting as indicated by Eq. (11). The distinction be- 
tween using Eqs. (11) and (12) for determining a climb 
technique is demonstrated graphically in Fig. 3. This 
figure shows contours of constant excess power per 
pound of weight superimposed on contours of constant 
specific energy. The two different paths shown in the 
figure result depending upon whether Eq. (11) or (12) 
isused as the governing condition. The use of Eq. (12) 
corresponds to the choosing of a speed where the excess 
power contours are tangent to a line of constant alti- 
tude. On the other hand, the use of Eq. (11) corre- 
sponds to the choosing of a speed where the contours 
are tangent to a line of constant specific energy. 

The advantages of the energy method in a minimum 


time problem are much more pronounced in the case of 
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a supersonic fighter, such as described in Fig. 2. Here 
the conventional method defined by Eq. (12) 1s vir 
tually useless because of the large speed change in 
volved if a speed near the maximum is desired at the 
The energy method, on the other 
The path defined by 


end of the climb. 
hand, is well adapted to this case. 
Eq. (11), which consists of selecting at each level of 
specific energy the speed and altitude for maximum ex 
cess power per pound of aircraft, is shown in Fig. 4. 
The flexibility of the energy method is here demon 
strated in that it may describe a climb with acceler 
ation or even an acceleration in a dive as part of the 
optimum climb path. The solution to this problem, 
essentially intractable by 
thus easily obtained by the application of the specific 


conventional methods, 1s 


energy concept. 

The effect of the difference in path on the time to 
the final speed and altitude can be visualized qualita- 
tively in Fig. 5. The time represented by the integral 
(S) is shown in this figure as the area under the curves 
of the reciprocal of the excess power per pound of 
weight between the limits of integration (/, lV’); and 
(#/W)s. These curves are plotted at a constant alti- 
tude so that the variation in specific energy implies 


also a variation in velocity along the curves. The ap- 
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plication of Eq. (11), which states that the minimum 
value of the integrand be selected at each abscissa 
(E/W), defines the envelope of the series of curves. 
The area under the envelope between (//W); and (£ 
W)s then represents the minimum value of the integral 
and, hence, the minimum time. The curve correspond- 
ing to the application of condition (12) is also shown 
and is seen to give a somewhat greater area. This 
comparison is not strictly valid because, in general, 
the final condition of speed and altitude desired will 
not lie on the envelope—for example, point B in Fig. 
5. The specific energy method gives the time from A 
to the point on the envelope that has the same specific 
energy as B—that is, point C. It essentially assumes 
that the zoom—that is, the conversion of the excess 
kinetic energy into potential energy necessary in going 
from point C to point B—can be accomplished in zero 
time. This, of course, is impossible, so that, in gen- 
eral, it will be necessary to depart from the envelope 
before reaching point C. The energy method generally 
results in the technique of acquiring a relatively high 
kinetic energy at low altitudes, where the excess power 
is the greatest for conversion, if necessary, into poten- 
tial energy at the higher altitudes. The modification 
of the integral due to departure from the envelope will 
depend upon the maneuvering capabilities of the air- 
craft and the amount of energy gained in the zoom. 
Some judgment will be required in this modification, 
although calculations indicate that its effect on the time 
in a long climb is small so that the greater portion of the 
theoretical advantage can be attained. Thus, the time 
from A to B via C (approximately), which is the mini- 
mum time path defined by the energy method, will be 
less than that from A to B along the conventional 
climb path. 

The problem of minimizing the fuel to achieve a speed 
and altitude can be treated by the energy method in the 
same fashion as the minimum time problem. This 
problem is significant in the case of aircraft, such as 
research vehicles, for which the flight durations at the 


desired test conditions are to be maximized. The fuel 
integral can be written as 
*(E/W)s | E 
W, = —— : a( : (13) 
ayw), 4@(E/W)/dW, \W 5 


But 


d(k,W) 
dW, 


d(k/W)/dt  (T — D)V 
dW,/dt  WTS.F.C. 


The only new variable introduced is the specific fuel 
consumption, S.F.C., which is a function of speed, 
altitude, and throttle setting. However, in a climb or 
accelerating condition, the throttle is generally wide 
open or at least at a fixed setting, so that the S.F.C. 
is a function only of speed and altitude. This permits 
us to write 
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ie m) ees 4 : 
dW, Ee oa (=) a 


The integral can now be written as 
(E/W): E E 
W, = f ( a v) a( :) or 
= arm, Mi i 
*(E/W): E E 

W, = fs ( i) a ( ) (15) 

j J, =a J5 W W ) 

The application of the calculus of variations yields the 


following condition for the path of minimum fuel ex- 
penditure: 


ra) eas 0 
oh WT S.F.C. E/W =const % = 


2 ay _ p)yv 
“ =0 (16) 
oV LWTS.F.C. \e/wacone ” 


The use of this result in an actual problem is completely 
analogous to that in the minimum time problem. . 
Since the primary purpose of the application of the 
specific energy concept to two typical problems was 
the demonstration of its usability, several details in the 
process were neglected. A-nong these are the variation 
of gross weight in the climb and the difference in the 
drag from the level flight value in the maneuvers re- 
quired to follow the optimum path. However, these 
details should offer no greater difficulty than that en- 
countered in the usual analysis of climb performance. 
Therefore, on the basis of this demonstration of its 
adaptability to problems for which the conventional 
methods yield no simple solutions, it is felt that the 
concept of specific energy deserves wider use in modern 
aircraft performance work. 
BALANCE APPROACH TO RANGE 
EQUATION 


THE ENERGY THE 


Range in the classic textbook approach is treated as 
a special problem completely independently of the other 
equations of performance. Here, however, the deri- 
vation of the range equation is made using the same 
energy balance relationship that was used to obtain 
This 


relationship as the starting point in the derivation is 


the rate of change of specific energy equation. 
rewritten here for convenience, 
Wd( kW) = mlldWw, — DV dt (2) 


since |’ dt = dR, Eq. (2) mav be written as 


weld... W /(E 
dW, — d 
D D W 
Making use of the relations for normal airplane cruising 
flight, 


dR = 
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dW = —dW, 
L =W 


, we may write 


(L/D) 


L\ dW i E 
dR = ~ null. ( )' = -( )a( ) (17) 
D/ "i D VW 


Now making the assumption that m and (LD) are 
constant during the cruising portion of the flight, Eq. 
(17) may be integrated to give 


P ‘i (5) W, (EVF(2) (X)] " 
_ oll; —_ _ (IS) 
” Dp)" w, D/\\w/, w/, 


Recalling that m//, represents the useful work of pro- 
pulsion done per pound of fuel, it may be expressed 


so that D = 


with the use of Eq. (4) as 
rv rv V 
nll. = = ~ Dana Aa 
(dW,/dt) TS.F.C. S.F.C. 
Substituting this relation into the first term of Eq. (18) 
and discarding the second term, the more familiar ex- 
pression for range is obtained 


R = (V/S.F.C.) (L/D) In (W,/W;) (19) 


This equation is recognized as the usual expression for 
the range of a jet-powered aircraft. The features dis- 
tinguishing Eq. (18) from this more familiar form of 
the range equation are its second term, as well as the 
parameters appearing in the first term. 

It may be seen from Eq. (18) that the need for the 
additional term arises only if there is a change in energy 
level during the cruising portion of the flight. Since 
the constant (L/D) operation, assumed in order to per- 
form the integration, essentially implies a constant air- 
plane attitude, the dynamic pressure will decrease with 
the gross weight. This type of operation will therefore 
generally result in a change in energy level and, strictly, 
will require the consideration of the additional term if 
the cruise is considered apart from the other portions 
of the flight. 

The magnitude of this term may be readily ascer- 
tained for the case of a turbojet-powered aircraft. 
For operation under the conditions of constant throttle 
setting, constant Mach Number, and altitudes above 
35,000 ft., the following relationship between a change 
in altitude and a change in gross weight may be written 


(see Appendix) : 


dW/W = —(g/RoIo)dh (20) 


Since a constant Mach Number at altitudes above 
35,000 ft. implies a constant speed in a standard at- 
mosphere, the energy level changes only by virtue of an 
altitude change, so that Eq. (17) may be written as 


ste ical < () ‘ : 
( = —m D W D an (21) 
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Using Eq. (20), 


L\dW g /L\aw 
1iR = aan iT, ) - ( ) 
: . & Ww 7 RT \D) W 


Combining 


g L\ dW 
dR = —\ eH. - => ( 
RoTv.1 \D/ 


and integrating, using the same assumptions as before, 


g | (=) W, 
rs a (22) 
RoTo_}) \D W, 
V J 8 Wo 
r=| a | ( ) 1m ( ‘) (9° 
SFC. RoT D W, 


In this form of the range equation, the second term ap- 
pears as a correction on the J’/S.F.C. factor. 
turbojet having a S.F.C. of approximately 1.0 and 


R= Ez 


or 


For a 


cruising at a speed of 500 m.p.h. at altitudes above 
35,000 ft., the correction on the range is approximately 
0.8 per cent. Although this analysis shows it to be 
small, a correction for the climb during cruise is ap- 
plicable. The use of Eq. (18) or the more specialized 
form, Eq. (22), should facilitate the making of this cor- 
rection. 

Where the additional term can be neglected, Eq. (18) 


may be written as 


R = mll, (L, D) In (We Wy) 


(24) 


The unique feature of the range equation in this form 
is the isolation of the fuel characteristics from the char 
acteristics of the power plant. 

The heat content of the fuel, //,, is a property of the 
fuel only, whereas the specific fuel consumption used in 
the more familiar forms of the range equation is de- 
pendent upon the characteristics of the power plant, as 
well as those of the fuel. For hydrocarbon fuels such 
as gasoline, this heat content has a value of approxi 
mately 18,800 B.t.u.’s per Ib. 
the B.t.u.’s to foot-pounds of energy using the mechani 
cal equivalent of heat, it can be seen that the dimension 
is that of a length, so that /7, can be expressed in the 
This value for 


However, by converting 


units in which range is usually given. 
gasoline is about 2,760 statute-miles. The dimension 
and the order of magnitude of aircraft range are there 
fore given by the heat content of the fuel. This fact 
makes clear the direction in which large improvements 
in range can be made. The advantage of nuclear fuels 
lies in their enormously high heat content. On the 
other hand, the poor range performance of rocket air 
craft is largely the result of the low energy content of 
itsfuel. ‘‘Fuel,’’ as used here, is considered to be made 
up of the oxidizer a rocket aircraft necessarily carries, 
as well as its actual fuel. The heat content of the fuel 
and oxidizer combination is as low as one-fifth of that 
of fuels such as gasoline used in air-breathing engines. 




















194 JOURNAL OF THE 
OVER-ALL 
POWERPLANT EFFICIENCY 
5 T T 
ALT = 30,000 FT 
"| COMPOUND TURBO 
3| ——_S _ TuRBOvET__— \RAMJET 
TURBOPRO | ili 





\\ 


| 7 RAMJET 
d a 


ie) 5 1.0 1.5 2.0 2.5 
MACH NUMBER 


Overall power-plant efficiency comparison. 


























Fic. 6 


The other parameter of Eq. (24) not found in the 
more familiar forms of the range equation is the over- 
all power-plant efficiency, m. Considered in the light 
of this parameter, the power plant is simply a device 
for converting the chemical energy contained in the fuel 
into useful work of propulsion on the aircraft. The 
efficiency with which it accomplishes this conversion is 
the overall power-plant efficiency. A typical value 
for this efficiency for, say, a reciprocating-engine- 
propeller combination operating near its peak, is 30 
per cent. As shown before, the product of this effi- 
ciency and the heat content of the fuel gives the useful 
work done per pound of fuel. 

The form of the range equation given in Eq. (24) is 
not particularly convenient for purposes of range calcu- 
lations of specific airplanes, since the power-plant char- 
acteristics are usually given in terms of specific fuel con- 
sumption. The advantage of the alternate form lies 
rather in its clear isolation of the four basic elements 
that determine aircraft range. For this reason, as well 
as for its general applicability to all power-plant types, 
this form is well adapted to range problems of a general 


nature. 


RELATIVE RANGE CAPABILITIES OF AIRCRAFT AT 
SUBSONIC AND SUPERSONIC SPEEDS 


Eq. (24) is particularly appropriate for the investiga- 
tion of the relative range capabilities of aircraft at sub- 
sonic and supersonic speeds. This problem has been 
considered elsewhere in detail, so that the discussion 
here will be limited simply to an interpretation of the 
results in terms of the four basic elements of aircraft 
range. 

(1) //., the heat content of the fuel, is a function of 
the fuel only. Assuming that the same fuel is used at 
subsonic and supersonic speeds, this factor will intro- 
duce no difference in the range capabilities in the two 
speed régimes. 

(2) m, the overall power-plant efficiency, is a func- 


A plot 


tion primarily of speed and power-plant type. 
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of this factor in Fig. 6, calculated from data given in 
reference 3, indicates that essentially the same effi- 
ciency can be achieved at supersonic speeds as at sub- 
sonic speeds, assuming the proper choice of power plant. 
This factor need not, therefore, introduce any differ- 
ence in the range capabilities in the two speed régimes. 
An interesting point here is that, since I’/S.F.C. = 
noll,, it is to be expected that the specific fuel consump- 
tion will increase with speed for any one fuel unless 
engines become basically more efficient. 

(3) L/D, the aerodynamic lift to drag ratio, is low 
at supersonic speeds because of the fundamentally 
different character of the airflow. In fact, its maxi- 
mum value now appears to be approximately only one- 
fourth of that attainable at subsonic speeds. The re- 
duction in lift to drag ratio will therefore introduce a 
corresponding reduction in the range capabilities of 
supersonic aircraft. 

(4) Wo/W,, the initial to final gross weight, is a 
function of the fuel weight to gross weight ratio of the 
airplane. In turn, this ratio depends upon the frac- 
tions of the gross weight required for power plant and 
structure. It appears, at present, that these fractions, 
in general, are not far different at supersonic than at 
subsonic speeds for airplanes of a similar function. 
This observation is made in the case of the structure 
weight on the basis that the weight increasing effects 
of aerodynamic heating, higher fuselage fineness ratios, 
and higher wing span to root thickness ratios are offset 
by the reducing effects of the higher wing loadings. 
Similarly, in the case of the power plant, the effect 
of the lower (L/D) on the power-plant requirement is 
approximately offset by the lower weight to thrust ratios 
available in the jet engines suitable for propulsion at 
supersonic speeds. A reasonable tentative conclusion 
would thus say that the weight ratios at supersonic 
speeds, in general, will be similar to those at subsonic 
speeds. No large difference is expected in the range 
capabilities in the two speed régimes as a result of the 
weight ratio. 

The overall conclusion from the foregoing exami- 
nation of the parameters of the range equation is that 
the range capabilities in the two speed régimes differ 
primarily only as the result of a difference in the lift to 
drag ratios. Since this lift to drag ratio at supersonic 
speeds at the present time reaches a value only one- 
fourth of that attainable at subsonic speeds, the range 
capabilities are lower by essentially the same factor. 

The form of the range equation, Eq. (24), has per- 
mitted a clear interpretation of the results of the prob- 
lem discussed above. The use of this equation has 
thus been shown to be appropriate for range studies of a 


general nature. 


CONCLUDING REMARKS 


The energy approach to the general aircraft per- 
formance problem has yielded several important re- 
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ENERGY APPROACH TO 
sults. In climb performance, it has pointed to the use 
of specific energy rather than altitude as the signifi- 
cant independent variable. This concept has been 
shown to be particularly useful in solving minimum 
time and minimum fuel problems of high-speed air- 
craft. In range performance, it has provided a clear 
isolation of the elements that determine aircraft range. 
In the case of both climb and range, it has led to equa- 
tions incorporating corrections that have been found 
necessary in modern performance work. Lastly, it 
has yielded a basic energy balance relationship that 
has general application in aircraft performance anal- 


ysis. 
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Appendix 


ANALYTICAL EXPRESSION RELATING (CHANGE _IN 
FLIGHT ALTITUDE TO CHANGE IN GROSS WEIGHT 


A turbojet-powered aircraft operating under the 
common conditions of constant attitude and altitudes 
above 35,000 ft. maintains as well, approximately, the 
following conditions: (1) constant throttle setting and 


(2) constant Mach Number. The altitude in this type 
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of operation then depends only on the gross weight in a 
manner that can be analytically expressed as shown 
below. 

In level unaccelerated flight, 


L = W = (y/2)pM*SC, Al 
When C, and / are constant, a differentiation yields 
dW = (y/2)M*SC, dp 
or 
dW/W = dp/p (A2) 
Now, combining the atmospheric pressure law, 
dp = —pgdh (A3) 
with the perfect gas law, 
pb = pRoTs (A4) 
the following relationship is obtained: 
dp/p = —(g/RoTo)dh (A5) 
Equating (2) and (5), 
dW/W = —(g/RoTy)dh (A6) 
Above 35,000 ft., the temperature is constant so that 
dW,/W = —dh 20,940 (A7 


This expression then gives the analytical relationship 
between a change in gross weight and a change in flight 
altitude for the special conditions given above. In 
integrated form, this relationship may be written as 


ho — h, = 20,940 In (W, W. 











The Bursting Speed of a Rotating Plastic Disc’ 
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ABSTRACT 


The paper presents an analysis of the stresses and strains in a 
fully plastic, rotating, annular disc that has initially uniform 
thickness and is made of a strain-hardening material. This anal- 
ysis is based on Tresca’s yield condition and the associated 
flow rule and assumes that the elastic strains may be neglected 
in comparison with the finite plastic strains that are considered. 
The bursting speed of the disc is expressed in the form of a defi- 
nite integral that involves the strain-hardening function of the 
material. In general, this integral must be evaluated numeri- 
cally, but analytical evaluation is possible for certain strain- 
hardening functions. In particular, it is shown that, for linear 
strain-hardening, instability can occur only at the onset of plastic 
flow, whereas for logarithmic strain-hardening considerable 
plastic deformation of a stable character may occur before the 


process of deformation becomes unstable at the bursting speed 


NOTATIONS 


ay, by = initial inner and outer radii of disc 

er, €@, @: = instantaneous (true) rates of strain in the radial, 
circumferential, and axial directions 

hy = initial uniform thickness of disc 

h = thickness of disc during plastic flow 

k,/,m = constants in logarithmic stress-strain relation, Eq. 
(16) 

ry = initial radial distance of a particle 

r = radial distance of a particle initially at 7 

Sr, 8g, Ss: = nominal stresses in the radial, circumferential, and 
axial directions 

‘7 = critical nominal stress (initial yield stress ) 

u = radial displacement of a particle 

v = radial velocity of a particle 

a => by ao 

€r, €&, €: = nominal strains in the radial, circumferential, and 
axial directions 

«* = critical nominal strain 

n = u/ao 

g = 1ro/do 

p = density of disc material 

or, 09, 07 = true stress in the radial, circumferential, and axial 
directions 

o = twice the critical value of the true shearing stress 

w = angular speed of the rotating disc 


(1) INTRODUCTION 


— ANALYSIS OF STRESSES AND STRAINS in elastic 
or partially plastic rotating discs has been dis- 
cussed repeatedly.'~* In these problems the strains 
are of the order of magnitude of elastic strains and may 
therefore be treated as infinitesimal. As a _ conse- 
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t Assistant Professor of Applied Mathematics 

** Professor of Applied Mechanics. 


quence of this, the equations of equilibrium may be 
satisfied for the undeformed rather than the deformed 
Considerable simplification of the mathematical 
It would be 


disc. 
work results from this approximation. 
unrealistic, however, to treat the strains as infinitesimal 
when the bursting speed of a strain-hardening plastic 
disc is to be determined. In fact, as the angular speed 
of such a disc is gradually increased, considerable plas 
tic deformation of a stable character may occur before 
the process of deformation finally becomes unstable 
at the bursting speed. 

Since the prediction of the bursting speed of a strain- 
hardening plastic disc requires the consideration of 
finite plastic strains, the analyst who has to forego the 
simplifications resulting from infinitesimal _ strains 
will look for other possibilities of reducing the mathe- 
matical work. Thus, Zaid’ uses a deformation theory 
of plasticity rather than a flow theory, in spite of his 
realization of the questionable value of this type of 
theory. The alternative analysis offered in the pres 
ent paper has recourse to other simplifying assump- 
tions. First, it is assumed that the elastic strains are 
negligible in comparison with the plastic strains. 
Secondly, it is assumed that the material of the disc 
obeys Tresca’s yield condition’ and the associated flow 
rule. 

CONDITION AND FLOW RULE 


(II) Tresca’s YIELD 


In a rotationally symmetric stress field, the prin- 
cipal stresses are the hoop stress og, the radial stress 
o,, and the axial stress o,.. The principal shearing 
stresses therefore are (1/2)\o, — a,|, (1/2)\e, — a,), 
and (1/2)!o, — o,;.. According to Tresca’s yield con- 
dition, none of these principal shearing stresses can 
exceed a critical value that depends on the stage of 
strain-hardening of the material. Moreover, for plas- 
tic flow to occur, at least one of the principal shearing 
stresses must have the critical value. 

If only one of the principal shearing stresses has the 
critical value, the flow rule associated with Tresca’s 
yield condition stipulates that the instantaneous strain 
rate corresponds to pure shear in the plane of maximum 
shearing stress. The sense of this shear deformation 
must be i.e., it must correspond to 
the sense of the maximum shearing stress. 


“appropriate” 


If two principal shearing stresses attain the critical 
value, the flow rule admits any strain rate that can be 
considered as resulting from the superposition of ap- 
propriate states of pure shear in the two planes of 
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critical shearing stress. It should be noted that in both 
cases the volume is preserved during plastic flow. 

To express this yield condition and flow rule analyti- 
cally, denote by o 2 the critical value of the shearing 
stress for the considered state of strain-hardening, and 
by é, e,, and e. denote the instantaneous rates of ex- 
tension in the circumferential, radial, and axial direc- 
tions. For the state of generalized plane stress oc- 
curring in a thin rotating disc, the axial stress co, is 
zero. Moreover, for an annular dise of initially uni 
form thickness, it may be tentatively assumed that the 
hoop stress og and the radial stress o, are tensile and 
that o, everywhere exceeds ¢,. Under these circum- 
stances, Tresea’s yield condition and the associated 


flow rule lead to the following relations: 
fn = CS 8, > 
e-=0, & = —e,>0 (2) 


Since og and e, are thus known explicitly, the number of 
unknown functions is reduced and the analysis is 
thereby simplified. 

To complete the description of the assumed me- 
chanical behavior of the disc material, it is necessary 
to formulate a law of strain-hardening. For simple 
tension, let the stress-strain diagram have the general 
shape of the line OAB in Fig. 1. The material remains 
rigid until the tensile stress s reaches the value s*; dur- 
ing the ensuing plastic flow, the stress increases with 
the strain € according to 


[= f(e) (.5) 


The rate of strain-hardening ds/de has the initial value 
* e* (see Fig. 1) and decreases monotonically with 
increasing strain. 

In Eq. (3), s and ¢ should be interpreted as the con- 
ventional stress and strain computed with reference to 
the original dimensions of the test specimen. On the 
other hand, the quantities o», o,, 7,, o, €, €,, and e, 
introduced above should be interpreted as the true 


stresses and strain rates; they are defined with respect 


ROTATING 


PLASTIC pisc 10 


to the dimensions in the considered state of deforma 
tion. 

It remains to discuss the application of the strain 
hardening law, Eq. (3), to the states of stress and plas 
tic flow occurring in a rotating disc. The state of plas 
tic flow described by Eqs. (2) is independent of the value 
of the intermediate principal stress ¢ \loreover, 
this type of plastic flow is possible under Tresca’s flow 
rule when the state of stress is simple tension in the 
It is 
can be applied to the problem 


circumferential direction therefore reasonable 
to assume that Eq. (3) 
on hand, provided that s is interpreted as the con 
ventional hoop stress sy and ¢ as the conventional hoop 
strain €,. 


(III) Furry Piastic Dis« 


Consider a rotating annular dise of the uniform int 
tial thickness /», the initial interior radius a,, and the 
initial exterior radius )). From the elastic stress anal- 
ysis it is known that, as the angular speed w of the 
disc is gradually increased, the yield limit is first 
reached at the interior surface. At somewhat higher 
speeds there will be an inner plastic region surrounded 
by an elastic region. Any flow that may occur in the 
plastic region would have to satisfy Eqs. (2). If v 
v(v) is the distribution of the radial velocity at a generic 
instant during this plastic flow, the radial strain rate 
is €é, = Ov/Or. 
Eq. (2), the radial velocity must be independent of r. 


Since this must vanish according to 


Moreover, since all elastic strains are neglected, the 


radial displacement uw at the elastic-plastic interface 


is zero. This means that and, hence, the circum 


ferential strain rate eg, = vr must vanish throughout 


the plastic region. This region thus remains rigid 
until the elastic-plastic interface reaches the exterior 


surface of the disc. The angular speed w, at which this 


occurs is given by see reference 6, page 103, 


Eq. (9.58) | 
w)" = (238 p) 
The fully 
plastic flow 
9.59 


where p is the density of the dise material 
plastic stress distribution at the onset ol 
is readily analyzed {see reference 6, page 104, Eq. 
and Fig. 9.5]. It is found that, at this instant, the 
radial stress nowhere reaches the vield stress regardless 
of the values of b) > a) > 0. Thus, Eqs and (2 
apply, initially at least, during the ensuing plastic flow. 

Let 7) denote the initial radius of a particle that is 
found at the radius 7 when the speed is w. If the radial 
displacement is denoted by u, then 


r=>°Y7roT ” 


Here, r and u may be considered as functions of the 


independent variables 7) and w* It is convenient to 
use w? rather than w, since the sense of the plastic flow 
is independent of the sense of rotation or the sign of 


the angular speed. 
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As has been shown above, the first Eq. (2) requires 
that all particles have the same radial velocity at any 
This implies that the radial displace- 


given instant. 
Thus, the 


ment depends only on w? but not on 7p. 
material bounded in the undeformed state by the co- 
axial cylinders of radii 7) and ry + dry undergoes the 
radial displacement « and is then bounded by cylin- 
ders of radii r = r~% + u andr + dr = r + dry in the 
Simultaneously, the original thick- 
Because of the incompressi- 


deformed state. 
ness /ip decreases to /. 
bility of the material, hyrpdr, = hrdr or, since dr = 
dro, 


h = ho(ro/r) = hol[ro/(ro + u)] (6) 


The true radial stress o, is transmitted across a sec- 
tion that is proportional to hr. Since hr = hor, by 
Eq. (6), the conventional and true radial stresses have 
The true circumferential 


the same value: o, = 5,. 
stress o, is transmitted across a section that is propor- 
tional to fdr; the corresponding value in the unde 
formed state is Aydry = Iydr. The conventional cir 


cumferential stress is therefore given by 
Se = ao(h/ho) (7) 


This stress is related to the conventional circumfer- 


ential strain 


by an equation of the form of Eq. (3). 
. * 
The equation of equilibrium in the deformed state is 
O(rho,)/Or = hog — hpw°r* (9) 
With reference to the undeformed state, this equation 
may be written as follows: 


Sg — puwro(to + U) 
= /(uU 7%) 


O(/oS,) Ory 
pw"ro(ro + Ut) (10 
vanishes at the interior and exterior surfaces 
the integral 


Since s 
of the disc—i.e., for 7) = dp and 7% = J 
of the right-hand side of Eq. (10) between the limits 
With a = bo/ao, » = U/ad, 


and & = 79, do, this condition yields 


ad) and & must vanish. 


pay-w*/6 = / f(n/&) dé/[2(a® — 1) + 3n(a? — 1)] 
J 1 
(11) 


For a known strain-hardening function / and given 
initial dimensions, the right-hand side of Eq. (11) 
must be evaluated, analytically or numerically, for a 
set of values of 9 = w/a. Each such evaluation fur- 
nishes one point of the plot w* vs. vu. For 7 = 0, in 
particular, /(0) = s* and the integral in Eq. (11) equals 
(a — 1)s*. Eq. (11) therefore yields the correct value, 
Eq. (4), for the angular speed at which plastic flow 
begins. 

As long as an increase in w* is required to produce an 
increase in “, the considered plastic flow is stable. The 
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bursting speed w. corresponds to the maximum of the 
curve w” VS. u. 

The preceding analysis is based on the assumption 
that the radial stress o, is smaller than the circumfer- 
ential stress og. A formula for the maximum radial 
stress occurring at a given speed w can be obtained by 
carrying out the differentiation on the left-hand _ side 


of Eq. (10) and equating Os, 07, to zero. Thus, 


max. 5S, = Sg — pwro(%o + U (12 


where 7) is now the initial value of the radius at which 
the maximum of s, occurs. Eq. (12) shows that even 
this maximum of s, is smaller than the corresponding 
value of sy. Now, o, = s,and og > sy by Eq. (7). The 
maximum of o, is therefore smaller than the corre 


sponding value of ap. 


(IV) LINEAR STRAIN-HARDENING 


As a rule, the integral in Eq. (11) must be evaluated 
numerically, but strain-hardening laws can be devised 
that make analytical evaluation possible. The sim- 
plest law of this kind is represented by the line OAC of 
Fig. 1. With 

f(e) = s*[1 + (€/e*)] (15 
Eq. (11) yields 


pay*w* a — 1 + (n/e*) log a 
6s* 2(a*® — 1) + 3n(a? — 1) 
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Differentiation of Eq. (14) with respect to » shows that 
d(w*) dn = O when 


e* = 2(a? + a + 1) log a/3(a? — | (lo 


Eq. (15) establishes the critical rate of linear strain- 
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FIG. 5 








the value of Eq. (15) if plastic deformation is to_ be 


stable. 1S plotted 


In Fig. 2, e* as given by Eq. (15 


versus a. As a— |—1.e., as the width of the annulus 


* 


tends toward zero— the critical value of «* tends toward 
|. This result is familiar from Laszlo’s paper 

The sign of d(w*) dn as found from Eq. (14) does 
not change with increasing 7. This means that, for 
linear strain-hardening, instability can occur only at 
the onset of plastic flow. A disc that starts to flow in 
a stable manner will continue to do so as the angular 


speed is increased. 


(V) LOGARITHMIC STRAIN-HARDENING 


Another strain-hardening law for which the integral 


in Eq. (11) can be evaluated analytically is given by 


€ k log (1 + me 16 
These may be chosen 


where &, /, and m are constants. 


hardening. For a given value of a = by) ad, the strain- so as to fit an experimental stress strain curve With 
hardening parameter e* (see Fig. 1) must not exceed — the function (16), Eq. (11) yields 
pay "lw (al + mn) log (al + mn) — (1 + mn) log (/ + mn) — al log a . 
= ‘ 
OR Pia — 1) + 3n(a- — | 


Fig. 3 shows the actual stress-strain curvet for Al 


ss: 


+ The authors are indebted to Prof. E. D’'Appolonia, Carnegie 
Institute of Technology, Pittsburgh, for the use of this curve 


24S-T4 in simple tension (full line) and the curve ob- 
tained from Eq. (16), with k = 9,780 Ibs. per sq.in., 


= 77, and m = 8,060. Fig. 4 shows the left-hand 
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side of Eq. (17) vs. 7 = u/do for a = bo/ay = 2. The 


bursting speed w. is obtained from the maximum point 
of this curve. As can be seen, the disc bursts for » = 
no = 0.25. Since » = u/d and ¢€ = u/r, € = 700/Po; 
since the maximum nominal hoop strain occurs at 7 = 
a, the disc would burst at a hoop strain of 25 per cent. 
From Fig. 5, it can be seen, however, that this hoop 
strain is beyond the range of validity of the expert 
mental stress-strain curve. A more realistic picture 
is obtained by integrating the right-hand side of Eq. 
(11) numerically, using the experimental stress-strain 
curve in Fig. 3 (full line) as s = f(e). 
* (where s* 


Fig. 5 shows a 
plot of pay*w*/bs is as indicated in Fig. 3) 
vs. n, for a = 2. Again, the bursting speed, w., is ob- 
tained from the maximurh point of this curve, which 
occurs at n = nm = 0.14, and thus at a maximum nomi- 
nal hoop strain of 14 per cent. This maximum strain 
occurs before necking takes place in the test specimen 
and thus is in the range of validity of the experimental 
stress-strain curve. The results just obtained illustrate 
the fact that, for logarithmic strain-hardening as well 
as actual strain-hardening of a similar type, there is a 
considerable range of stable plastic deformation before 


bursting occurs. 
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on 
The Calculation of Wall Shearing Stress from practically constant It is believed, however (see reference 4 
uy Heat-Transfer Measurementsin Compressible that such a condition could also be assumed for the compressible 
377 Flows case. As a consequence it can be shown that ¥,, = AK = con 
stant and Wann = Yor = Yr = O, so that the momentum equa 
Nick S. Diaconis tion reduces to 
‘ Lewis Flight Propulsion Laboratory, NACA, Cleveland, Ohio v -0 (3 
. December 9, 1953 p= ~— " ” 
The solution of Eq. (3) with the boundary condition y,(x, 0) = 
, ‘ — , 0, is 
ides, I’ HAS BEEN SHOWN by Ludwieg! that the wall shearing stress 
799 of a laminar or turbulent boundary layer in an incompressible y = (An?/2) + constant (4 
vy can be determined from a heat-transfer measurement at the a ‘ ; = ; 
flow a ee : ” mae Che energy equation simplifies to the following 
surface The instrument used in that investigation was essen 
J : i . 7" >\T ee he 2 \2 ~ 
tially a small, locally insulated, heating element embedded in the vnl CU/PY)Tan = (9 1) A? (Wns (0 
1) Tr . > . E . x ao = x . . bs = P 7 - . ? 3 
a test surface The size of the instrument was restricted by the Substituting Eq. (4) into Eq. (5) and introducing the new vari 
condition that the thermal boundary layer generated by the — ables ¢and 7 yields 
) heating element be contained locally within the laminar sub P 
) : ‘ me , : nl ye — Ty (vy — 1).M,°A (6 
til layer. In the present analysis Ludwieg’s theory for such an in 
AUK ¥ P ° 
strument is extended to compressible flow over an insulated flat where 
plate With the same limitations on the design and operation — = 
“ P ‘ ; f- WE rr ss, 7»= VAPrn (7 
of the instrument as mentioned above, it can also be assumed 
for compressible laminar and turbulent boundary layers that Introducing the similarity variable a (see reference 1), where 
only the flow in the immediate vicinity of the wall or the laminar = n/(9 (2 
° P . a >= we 
sublayer will be affected in the region of the heated element 
This assumption then permits the use of the laminar boundary results in the transformation of Eq. (6) into the ordinary differ 
layer equations as the governing equations for this. analysis for ential equation 
_ both laminar and turbulent boundary layers (d?T/da?) + 38a2°%(dT/da) = (y — 1) ,2A (9 
If the specific heat and Prandtl Number are taken to be con = - Dek 
aa : : . . . Since Eq (9) is linear, the temperature 7 can be expressed as 
stant and if a linear viscosity-temperature relation is assumed, 
the compressible equations become in dimensionless notation [=fT7,+ 7 (10 
and in a transformed coordinate system similar to reference 2 an 4 —— ; , 
— | transformed coordinate system s si where 7) is the temperature distribution when the plate is at the 
Vome ntuni: i . . < 
adiabatic wall temperature and 7% is the perturbation caused by 
Wn Van¥rc = y¥ (1 slightly heating a small segment of the wali (i.e., at the heat 
i transfer instrument) from 7, to 7 s3oth 7) and 7» satisfy 
4 ¢ t : ‘ - “ u 
Eq. (9), which can then be solved for the following boundary 
Yul - orl (1/Pr)7 Ply > 241," te (2 conditions: 
where a= O(ie, 7 = 0 7) = Tua 
, , ‘ “ser . dT, /da = 0 { 
y L(V ew *l L* C = stream function : ws 11 
é eo : ; . T» © Ty — Tea 
/ lad BF = absolute static temperature ratio z 
> a= (1.e.,y > Oand ¢—~0) ] 0) 
ee ” 
dy Qe ° 
} = = transformed normal coordinate 
0 ul where 
v* 1* = longitudinal coordinate Tag = adiabatic wall temperature 
( a/b ) (PR / Pw viscosity-temperature propor 7, = wall temperature with heating clement operating 
tionality factor ail . a 
— The solutions for 7; and 7> are 
i a characteristic length 
limensional ity Pa * 
aati nie: atteaiied Ta Tr (+ — 1)M.2K, | E w { \ an] de (12 
= conditions ‘n the free strean JV JW 
(Indep ndent variables used as subscripts denote partial dif T= To = Fest be = DBE = KL) 
ferentiation i 
a 
If Eqs. (1) and (2) are restricted to the laminar sublayer, cer cx . i, . da 
yer, . — 
. ‘ a . . ilies e dy | dw + (Tata — Tw) ~ (15 
tain simplifications can be made. From the low-speed results JO J0 : 
° rs Z 7 ‘ a - ( , 
of Laufer,’ it appears that the shear in the laminar sublayer is Jo | 
201 











202 JOURNAL OF THE 


The addition of Eqs. (12) and (13) then gives the complete solu- 
tion of Eq. (9). Note that the constant A, does not appear in 
this solution and need not be evaluated. The constant A, how- 
ever, is related to the wall shearing stress 7,.* through Eq. (4) 
and can be expressed as follows, where the subscript w represents 
wall values 


V r5*/put*\?  [U_*L* 
a (14) 
U.* 7 


The Nusselt Number is then defined as follows 


K = (Yan)o = 


+) 


or” 
- dx* 
- h* i* 0 oy" J y*=0 
Nu = = — a 
ze Toa* — Tw* 
lu,*L* [' for 
dé 
Po 0 On n=0 - 
oo a (15) 
Tua ~ T. 
where 
|] «= VE Pr @*/L*) 
/* = length of heating element in measuring instrument 


Substituting the expression for the temperature distribution 
in Eq. (15) gives for the Nusselt Number in compressible flow in 
the physical plane 

Nu = 0.807(rw*pw*Pr)'/? (1*/pw*) (16) 


Eq. (16) shows that the compressible Nusselt Number is ob- 


* and #.* 


tained from the incompressible form by replacing p, 
by pw* and pw*, respectively. Although the compressible tem- 
perature distribution in the sublayer has an additional term (a 
function of stream Mach Number) that does not appear in the 
expression of reference 1, this term does not contribute to the 
heat transfer. Finally, Eq. (16) shows that, even in the com- 
pressible case with the use of an instrument of Ludwieg’s type, it 
is possible to determine local skin friction values (7,,*) for an 
insulated flat plate from measurements of the heat transfer at a 
heated element in the wall. 

Liepmann and Skinner® have obtained essentially the same re 


sult by means of a physical argument 
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On a Perturbation Theory Based on the Method 
of Characteristics 


CC. Lin 

Professor of Mathematics, Massachusetts Institute of Technology, 
Cambridge, Mass. 

December 1, 1953 


T IS KNOWN THAT a perturbation theory in the physical plane 
for supersonic flow and wave propagation is not free from 
difficulties. A method for improvement has been proposed by 
Lighthill and Whitham! for dealing with outgoing waves. On 
the other hand, the characteristic parameters are the natural 
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Fic. 1. 


Interaction of two simple waves. 


independent variables to use for hyperbolic differential equa 
tions, and it may be expected that, in a perturbation theory d¢ 
veloped in terms of the characteristic form of the differential 
equations, such difficulties will not occur. 

The method of characteristics is well known, but it does not 
seem to be extensively used in connection with the perturbation 
theory.” In fact, in this case, an additional convenience can be 
easily introduced. Since the characteristics are nearly straight 
in the physical plane, the mapping between the characteristic 
parameters and the physical coordinates can be made an ap 
proximately affine transformation, which is easy to visualize 
The complication that sometimes arises with the hodograph 
method is avoided. 

On the other hand, when the hodograph method is useful the 
characteristics are fixed in the hodograph plane. Thus, the 
present method would also embody the advantages of the hodo 
graph method. 

Briefly, the method works as follows. If (x, y) and 2(x, y 
are the solutions sought, we introduce the characteristic param 
eters @ and 6 and seek the solution in the form 


“w= wu (a, B) + a (a, 8B) +. 

v =v (a, B) + ew” (a, B) 4+ | 
x = x (a, B) + e'” (a, B) +... 

y= y' (a, B) + ey (a, 8B) +. 


Thus, the solution is given by a parametric representation. The 
initial approximations to x and y are expected to be linear func 
tions of a and p. 

Several examples have been worked out by the present writer 
and Miss Phyllis Fox,* using the method just explained. Con 
vergence proof is given in certain cases. The reflection of a 
Prandtl-Meyer expansion from the plane boundary will now be 
described to show how the method works. It is clear that the 
divergence of the characteristics is important here, and hence the 
perturbation theory in the physical plane does not work satis 
factorily. 

The problem is equivalent to that of the interaction of two 
expansion waves of equal strength, and we shall therefore con 
sider the general case of interaction of two simple waves. The 
interaction region is bounded by four characteristics forming a 
rectangle in the (a, 8) plane. 

Along the two characteristics OA and OC, the solution is 
known; thus, we have a characteristic initial value problem. The 
choice of the scales of a and 8 is still at our disposal. It is im 
portant, for convenience of application, that this choice be prop 
erly made to simplify the mapping between the (x, y) plane and 
the (a, 8) plane. This can be done simply by taking x = @ 
along OA and x = 6 along OC. The given physical conditions 
along OA and OC should be presented as functions of x. It has 
been found convenient to describe the state of the gas by using 
the Mach angle A and the angle 6 for the direction of flow. Thus, 
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ifo = f(x along OA and @ = eg(x) along OC, the characteristic 


yalue problem is formulated 


6 = ef(a), x =a 


6 = eg(8), x =8 


as follows: 
along OA} 
along oc} 


The equations for @ and A can be integrated in closed form, and 


we have 
6 = elf(a) 4 g(B)] | 
(3) 
¥(A ) ¥(Ao) = el f(a) (B)|$ 
where 
¥(A) = A — (w/2) + (1/p) tan (pw cot A) | (4 
uw? = (7 — 1)/(y + 1) \ 
The equations for x and y then give 
vr =(at+ B) + en/2d) [Bf(a) ag(p)| 
vy = Na B) + €/(n/2) [Bf( a) + ag(B)| . 
"a "8 (9) 
e (m/2) f(a) da 4 2(8) dB 
0 0 J 
where 
2 f , ’ 
m = sec? Ag {1 + [1/¥(Ao)] § 
‘ (\ = tan Ao) (6 
n = sec? Ay {1 — [1/p'(Ao)]} f 


The reflection of a simple 
obtained by taking g(6) = 
a = 8, where @ also vanishes. 


for the pressure distribution 
A detailed account of the 


wave from a solid boundary can be 
1(B) 
We have then 


The boundary y = 0 becomes 


7 2ef(a) 


= 2 e( A)af(a) 4 


= al 


along the solid boundary 
method and the discussion of the re 








sults will be published later. 
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Experimental Velocity Profiles for Adiabatic 
Supersonic Flow of Air in a Tube 


Joseph Kaye,* George A. Brown,t and John J. Dieckmannt 
Massachusetts Institute of Technology 
December 1, 1953 


a THE PAST SEVERAL YEARS, a program has been under 
way to measure recovery factors and heat-transfer coeffi- 
cients for supersonic flow of air in a round tube. It has been 
assumed in the interpretations of the data obtained that, for 
smooth entrance into a tube, the flow of air at supersonic speeds 
is attended by the gradual growth of a boundary layer in the 
direction of flow and by the filling of the tube cross section with 
this boundary layer. In the last few months, measurements 
of the velocity profiles of adiabatic supersonic flow of air in a 
round tube, '/2 in. in diameter, by means of small impact pres 
sure probes, have demonstrated that it is possible to permit a 
laminar boundary layer to develop and grow over a long length of 


tube flow. These data, establishing the existence of a laminar 
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boundary layer over a long length of tube flow, permit a more 
secure foundation to be established for the comparison of theo 
retical work with experimental results and also permit a better 
comparison of tube flow with supersonic flow over a flat plate 

These measurements of velocity profiles were performed as part 
of a research program sponsored by the Office of Naval Research 
on measurement of heat-transfer coefficients for supersonic flow 
of air. The details of the apparatus used and the results ob 
tained are to be found in previous reports.!~® 

Several small impact pressure probes were constructed with 
the severe limitation that the probes should result in a minimum 
interference with the supersonic flow in the tube. The '/2-in 
diameter tube is approximately 25 in. long, and the probe was 
The 


details of two such successful impact pressure probes are shown 


to be inserted about 2.75 in. from the downstream end 
in Fig. 1. The probes were mounted on a traversing rig, which 
was situated in the downstream stagnation tank of the apparatus 
This traversing rig enabled the position of the tip of the probe 
to be varied, in small increments, from the center of the tube to 
within about 0.004 in. from the wall of the tube. An electronic 
device was used to indicate when the probe tip was in actual phys 
ical contact with the tube wall; this device operated very suc 
cessfully. Furthermore, the design of the rig permitted a rota- 
tion of the probe by 90° to allow measurements of the asymmetry 
of the velocity profiles for two planes at right angles to each 
other. The actual vibration of the probe tip in the experi- 
ments, especially near the tube wall, was found to be less than 


0.003 in. in amplitude. 
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1.0 mo— The preliminary measurements of velocity profiles presented lees 
° . . . whe 
here were obtained with a stainless-steel nozzle designed to pro- Reyn 
° = " e) 
gor duce a supersonic flow at a Mach Number of 2.8 at the entrance 
3 plane of the tube, 25 in. long. For adiabatic flow the inlet dj 
° 0.6 ADIABATIC FLOW PROBE J, VERTICAL ameter Reynolds Number could be varied from 0.3 to 1.0 X 103 in wh 
7 TRAVERSE AT STATION 17 to produce a laminar boundary layer over the entire length of to cc 
{ 
" 0.4 RUN CODE Rene, 210° this tube. For an inlet diameter Reynolds Number of about ing th 
a P-17 fe) 1.00 2 X 10°, a shock was established beyond the entrance of the tube: te 
3 0.2 ye . : 7 the resulting flow downstream of this shock produced a subsonic both 
a turbulent boundary layer extending over the rest of the tube Since 
- ; : ) 
0 length. Velocity profiles are shown here for both the laminar xcell 
e | 
O Ol O02 03 04 05 O06 O07 08 O09 10 andthe turbulent boundary layers. this 2 
RADIAL POSITION, y/a Fig. 2 compares the profiles obtained with probes // and J in CON 
FIG. 3-EFFECT OF INLET DIAMETER REYNOLDS NUMBER The velocity ratio is the ratio of computed local velocity to com theor 
ON VELOCITY PROFILES puted centerline velocity; the radial position is the ratio of the the 
. e a’ 
distance of the center of the probe tip from the tube wall to the tatio 
radius of the tube. The agreement of the results obtained with Fig 
10 Pon —— the two probes is excellent. a tel 
Fig. 3 shows the effect of changing the diameter Reynolds a 
+ -. a DC 
09 Number at tube inlet on the velocity profiles near the exit of the watt 
tube. Increasing this Reynolds Number from 0.8 to 1.0 X 105 . 
0B increases the slope of the velocity profile in the laminar boundary 
layer, but the nearly linear characteristic slope is evident in 
07 each of the three profiles. ; 
THEORETICAL PREDICTION FOR ADIABATIC Fig. 4 compares the velocity profiles obtained at station 17 
> os FLOW, LAMINAR BOUNDARY LAYER ON (2.75 in. from tube exit) in the tube with those predicted by Chap- vhere 
> fli pice gs gp (6) man and Rubesin® for adiabatic laminar flow over a flat plate flow 
© os / and with the measured prefiles of Blue? for flow of a laminar Eq. (+ 
mS / " cut ° : 7 
= / EXPERIMENTAL CURVE FAIRED THROUGH boundary layer over a flat plate. The abscissa is the customary A lirst 
aS / DATA OF BLUE (7), LAMINAR BOUNDARY . lent 
> 142 LAYER ON FLAT PLATE one used in analyses of flat-plate flows—namely, sccaiat 
— / Re, = 10 to 0.39 x 10° : tube | 
S f M= 20 n (y L) V Rex (1 betwe 
WwW 
> as-/ ; ; : , ; plate 1 
FY oO RUN CODE Rexi7 x 10° Mei7 where L is the distance from the leading edge, and the length ’ ne 
: oe ee ‘ig 
P-17 O 47 238 Reynolds Number, Kez, ts given by a ) 
tur 
02 P-18 ® 3.4 2.29 Re = " — 
P-15 a) 1.9 2.15 f : a veloci 
Ol where #o and vo denote the free-stream velocity and kinematic layer | 
viscosity, respectively. For the case of tube flow, however, tube-t 
fe) - where the centerline velocity and fluid properties vary with the than 1 
Oo 2 4 6 8 10 12 14 16 18 20 . : fortuit 
distance downstream from tube entrance, a convenient though ortul 
RELATIVE DISTANCE FROM WALL, 7 ; eae : Th 
arbitrary definition of » is used here—-namely, . 
FIG.4 COMPARISON OF VELOCITY PROFILES FOR TUBE FLOW AND , in the 
FOR FLOW OVER A FLAT PLATE n (y/L) V Ren; 3 
1.0 | : | oO—-00— ea 
' Measur 
— - = 
- o-O-. ——— hanics 
- 
5 O8F-F —_ - 
-* Fart 
= vA EXPERIMENTAL CURVES eral Di 
. eo) F STALDER, RUBESIN, TENDELAND (8) tion, AS 
O 06 TURBULENT —-: —M= 2.4 (FAIRED) : i 
; oong, 
oS LAMINAR ==*-——M=2.4 (CURVES) Sube 
ea 74, pp 
i PLATE | THEORETICAL CURVES ee 
- OFF FLOW | TURBULENT ——— —— ste 
5 FRANKL, VOISHEL (9) * Kay 
7 te LAMINAR amen anmemnenans Supe 
BW o2 e CHAPMAN, RUBESIN (6) M=2.4 Heat 1 
° — of Iowa 
S TUBE RUN CODE Reoci x10 Meiz7 ‘Cha 
0 | | | | P-24 ®O ZANZo 0.861 ols 
i de Septem 
o «4 6 ¢ © @ 6 6 8 2 22 Blu 
§ Stal 
RELATIVE DISTANCE FROM WALL, y/@ fe Kae 
very F 
9 Fra 
FIG. 5 - COMPARISON OF VELOCITY PROFILES FOR TUBE FLOW AND FOR Flat Pl 
TM No 
FLOW OVER A FLAT PLATE 








ited 
pro- 
ince 
t di 

10) 
h of 
0ut 
ibe: 
onic 
ube 


inar 


om 

the 
the 
‘ith 


ids 
the 
105 
ary 


in 


ip- 
ite 


lar 





READERS’ 


where L is the distance from entrance plane of the tube and the 
Reynolds Number is arbitrarily chosen for convenience to be 


Ren; = taly/vea (4) 


in which the fluid properties are to be evaluated at tube entrance 
(to correspond to flat-plate flow), or at station 1, and by assum 
ing that an isentropic core exists for the flow at tube entrance 

In Fig. 4, the tube-flow profiles show excellent agreement with 
both experimental and theoretical profiles for flat-plate flow. 
Since the definition of 7 used in this comparison is arbitrary, this 
The effect of 


this arbitrary definition will be examined later in greater detail 


excellent agreement may be, in part, fortuitous 


in comparing the velocity profiles with the results based on a 
theoretical analysis of tube flow. In Fig. 4 the values of 1/¢1; are 
the Mach Numbers computed for the centerline velocities at 
station 17 

Fig. 5 shows another means of comparison of velocity profiles 
for tube flow and flat-plate flow for both laminar and turbulent 
boundary layers. The abscissa is the ratio of distance from the 
wall to the momentum thickness of the boundary layer, defined 


5 
pl ul “ 

0 l dy (5 
Pi i it 


for plate 
Alt hough 


by 


vhere subseript 1 refers to the free-stream conditions 
flow and to the centerline conditions for tube floy 
Eq. (5) is not strictly correct for tube flow, it is used here only as 
i first approximation. For the laminar boundary layers, ex 
cellent agreement is obtained between the velocity profiles for 
tube flow and plate flow, although the remarkable agreement 
between the tube-flow profile and the theoretical curve for a flat- 
plate flow may be coincidental 

Fig. 5 shows also that, for subsonic flow in the tube, where a 
turbulent profile should be obtained after a shock, the measured 
velocity profiles agree well with those for turbulent boundary 
layer on a flat plate. Again the agreement of the experimental 
tube-flow profile with the theoretical flat-plate profile is better 
than with the experimental flat-plate profile; this result may be 
fortuitou 

The assistance of Edward A. Sziklas and W. S. Wu in aiding 
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The Numerical Solution of Certain Nonlinear 
Supersonic Flow Problems on the IBM Card 
Programmed Calculator 


R. E. Bower and J. H. Higman 

Research Engineers, Grumman Aircraft Engineering Corporatior 
Bethpage, L.I., N.Y 

December 8, 1953 


— COMPLEXITY OF NONLINEAR SUPERSONIC FLOW PROBLEMS 
precludes an exact analytical solution except in certain rare 
instances. Generally, the powerful numerical procedures, such 
as the method of characteristics, are resorted to, but their appli 
cation is restricted because of the extensive amount of time 
which must be expended by better than average computing per 
sonnel. The recent increase in the availability of automatic 
computing equipment, however, is eliminating this obstacle, and 
procedures previously considered prohibitive in practice now 
take on new elegance. Clippinger and Gerber! have already 
treated the automatic computation of bodies of revolution at 
zero angle of attack on the ENIAC. This note is concerned with 
similar work performed on the IBM Card Programmed Calcu 
lator (C.P.C The C.PC. is an electronic digital computor 
which, in contrast to the ENIAC, is available to most computing 
groups 

The study is restricted to steady-state problems that are func 
tions of not more than two independent variables. Future work 
is contemplated, however, to consider the effects of a third vari- 
able. The numerical procedures are in the main developed from 
the theory of characteristics for hyperbolic partial differential 
equations. In the case of conical flow, however, Busemann’s 
radius of curvature solution in the hodograph plane? is used, and, 
of course, for two-dimensional oblique shock waves, the conven 
tional methods for solving a cubic are employed 

The characteristic relations are revised to avoid trigonometric 
functions. This not only results in a saving in machine time but 
also allows more versatility in the use of the C.P.C. installation 
at Grumman in that problems can be card-programmed for us« 
with the more general purpose type of calculator setup which is 
capable of performing the basic arithmetical operations of multi 
plication, division, addition, subtraction, and square root. The 


form of the equations used is 


dr uv +a Vv V2 — a? | 
== _ (la 
dx u2? — a? m 
du di a?’ 
e a 4 
dr d uv +a V V? — a? 
rt a \ | 2 a?* dS 
« inst 0 (lb 
, yVR dn 
and 
d 7 — a? 
= =} Da 
dx uw+avi - a? 
du di a? : 
rr 4 = 4 
dx dx uv +aV V? — a? 
v aV V? —a’?dsS : 
€ + : = () Pb 
r yVR dn 


where u and v are the velocity components of V along the cylin 
drical or rectangular axes x and ¢, a is the local sound velocity, 
dS/dn is the gradient of entropy normal to a stream line, y is the 
ratio of specific heats, R is the universal gas constant, and e is 
zero for two-dimensional flow and unity for axially symmetric 
phenomena. The arrangement of Eqs. (1) and (2) avoids inde 
terminate forms when u = a. 

The solution of a particular problem by the method of char 


acteristics involves the repetitious application of certain basic 
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TABLE 1 

C.P.C. Running Times of Some Representative Calculation 
Decks 


-Running Time (Min.)— 
First* Each 
iteration 


Function of Deck approximation 


Calculation of an interior point in a rotational, 


axially symmetric flow field 3 0.6 
Calculation of the axially symmetric flow con 

ditions at a fixed point on a physical bound- 

ary (at a wall) 1.4 0.5 
Calculation of flow conditions at a point on the 

axis of symmetry 1.0 0.2 
Calculation of a point on a shock boundary in 

uniform, axially symmetric, rotational flow 2.1 0.9 
Calculation of a point on a shock boundary in 

a nonuniform, axially symmetric, rotational 

flow field 3.9 Me 
Calculation of the flow crossing an oblique 

shock wave 1.7t 


* Includes punching of results. 
+t Total calculation time. 


computations based on Eqs. (1) and (2) and appropriate bound- 
ary conditions. The programming on an automatic machine 
consists then in setting up these elementary calculation decks 
with a provision for their combination. The decks in this study 
were card-programmed for use with the C.P.C. storage capacity 
set at 48 positions. With such a limitation it is not possible to 
store sizable portions of the flow field so the effective storage of 
the calculator must be increased by the development of a system- 
atic method of punching and storing the space coordinates and 
velocities identified with each mesh point. As these quanti- 
ties are calculated for each point, they are punched with proper 
coding for future use and, when necessary, held in storage for 
immediate use in the solution of a subsequent mesh point. 

A total of 30 basic computation decks were programmed on 
the C.P.C. They allow the completely automatic solution of 
such nonlinear supersonic problems as (1) the flow transition 
through oblique shock waves, (2) conical flow phenomena for 
circular cones, (3) the rotational flow over two-dimensional 
bodies, (4) external irrotational or rotational flow about pointed 
or open-nose bodies of revolution, (5) the irrotational flow through 
diffusers or nozzles, and (6) jet exhaust flow. The programming 
of an extremely lengthy computation such as an axially sym- 
metric shock boundary in a nonuniform flow field was attempted 
as perhaps a practical upper limit for a single calculation per 
formed on the C.P.C. Although the solution was successful, 
the complexity and tediousness of the programming indicated 
the impracticability of treating a more involved problem with a 
single deck. The development of the nonuniform shock bound- 
ary computation, however, extended the list of soluble flow prob 
lems to include ducted two-dimensional or axially symmetric 
bodies. 

Because of the obvious enormity of the problem, a detailed 
error study is not considered) The characteristic relations are 
solved by a second-order iterative type of solution, and eight 
digits are carried in the fixed decimal system to minimize round 
off errors. The error considerations of reference 1 allow in- 
creased accuracy for a given amount of time in some cases by an 
extrapolation to zero mesh size. With the assumption that the 
error varies as h? = 1/n?, where n is the number of integration 
steps taken on the body or along an initial characteristic, Rich 
ardson’s h? extrapolation coefficients may be used on two solu- 
tions performed at different grid sizes. 

The saving in time realized by the C.P.C. is dependent upon 
the complexity of the problem and the type of hand method em- 
ployed in the comparison. Perhaps quoting some actual ma- 
chine running times for the calculations as is done in Table 1 is 
The number of iterations necessary for con 


the best measure. 
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vergence depends upon the Mach Number and grid size, byt 
usually two or three are sufficient. The time to solve any par- 
ticular problem is easily obtained from the number of mesh 
points to be considered, since the solution is automatic and con- 


tinuous. 
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Effect of Local Modifications in Redundant 
Structures 


E. Kosko 
Assistant Chief Stress Engineer, A. V. Roe Canada Limited 


December 9, 1953 


oo APPEARS TO BE A NEED for some means of eval 
uating the effect of a local modification in a redundant struc 
ture on the stress distribution and on the deflections. Refer 
ence | suggests that an extension of Best’s ‘‘Stress Area Method’? 
is suitable for that purpose. A more recent contribution® points 
out that a similar condition exists in electrical engineering, where 
the so-called Compensation Theorem provides a solution. This 
writer, however, being unfamiliar with electrical networks, has 
found the description of the method somewhat difficult to follow 
It is therefore proposed to present an alternative formula with a 
derivation based on matrix algebra. ' 

Let us assume that a linear relation exists in the form 


A,oU =P (1 


between 2” load components 1P,, Po, ..., Pa} forming a column 
vector P and the corresponding deflections }| U1, Us, _ Uf ofan 
elastic system, where Ay is a square matrix of order n, generally 
symmetrical, called the elastic or stiffness matrix. For any sys 
tem encountered in engineering, this matrix is certainly non 
singular, and we may assume that a unique solution of Eq. (1) 
has been obtained in the form 

U = A, 'P 2) 
The elements of the inverse matrix Ay | are the influence coeffi 
cients of the structure, and 49! is called the influence or flexibility 
matrix. 

Assume now, as in reference 3, that the modification in the 
structure affects a square submatrix of A) ~! which is of the order 
m <n. Without loss of generality, this submatrix may be placed 
in the top left-hand corner of the stiffness matrix; such a position 
can always be obtained by suitable interchange of rows and col 
The modified stiffness matrix may then be written 


A = Apo +a (3 


b O 
= (4 
“ [f 4 


b being a square matrix of order m, while the symbols O represent 
zero submatrices. In order to invert the modified matrix A, it 
is advantageous to present it in the form of a product rather than 


umns 


where 


a sum—viz., 


A = QAy 5 
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where 

Q =I, + aAyg™ (6) 
Then it is possible to write 

A>! = Ay-'Q-1 (7) 
The ‘‘quotient’”’ matrix Q is of a special type and is generally 
easier to invert than A. Let the leading m by m submatrix of 
4,)~! be denoted by K, so that Ay~! is partitioned as follows: 


Pre a . 
aie 2 9) 


where L, 1/, and N are (m by r), (r by m), and (r by r) sub 


matrices, respectively, withr = —m. Then 


bO KL bK bL 
O0= f. + 4 = 7. + = 
OO MN OO 
[i+ bK a os 
O I, 


where /,, and /, are unit matrices of order m and r, respectively 


Using the notation 
g=In+ bk (10) 


the inverse of Q is obtained 


qt —q oL 
1 ) 
QuU= | O 1, | (11 


provided g is nonsingular. It is thus seen that, instead of in 
verting a matrix of order 1, it is necessary only to invert a matrix 
) of lower order m, using the results of a previous inversion of the 
original matrix of order x. Substituting Eqs. (8) and (11) into 


Eq. (7), we obtain the desired result 


Kq 1 (Im — Kq a (12 
Mq7', N — Mq dL ae 


A simplification results if we consider the difference between the 
inverse of the original matrix and that of the modified one. It 
iseasy to verify that this correction term may be put in the form 


ofa matrix product 


A=Ay!-—-At= ‘le (KL) (13) 


Had a ‘“‘quotient’’ matrix been introduced as in Eq. (5), but 
on the right-hand side of Ao, as follows: 


A = AoR (5’) 


with 


the result would have been obtained in a slightly different form 


i. KL} 3’) 
ro br “[|KL} (13 


y=/fIf,, + Kb (10° 


namely, 


where 


It is to be noted that, in case of both A and > being symmetrical 
ris the transpose of q. 

When only a single term of Ao is modified, assuming this to 
be the leading term of the principal diagonal, we have 


gq = 1+ dk 
i scalar), where & is the leading term of A,~', and 
q 1h = b/(1 + Dk) 


lhe correction matrix becomes now 


b k 
A= [ | [k, L] (14 
1+ bk LM 


L| is the first row and |k, A/} the first column of Ay! 


where [fk 
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This correction is seen to be a matrix of rank 1, each row of which 
is proportional to the first row and each column proportional to 
the first column of Ay~'. In the more general case of Eq. (13), the 
correction matrix is also singular, but its rank may be as high as 
m. 

The formulas derived above may also be used when only an 
approximate answer is required. Noting that the elements of 
the product matrix bX will generally be small compared to unity, 
the inverse of g may be expanded into a matrix power series 


g-' = In — bK + (bK)? — (bK)§ + 


which in many cases should converge rapidly. The series may 
be broken off at a convenient number of terms, and their sum may 
be used instead of the exact value of g~!. 

As with the Compensation Theorem, the above formulas are 
applicable to matrices of any kind. It has been thought expe 
dient, however, to limit the discussion to structural matrices, 


which today are the subject of wide interest. 
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On the Calculation of Laminar Boundary-Layer 
Stability 


Martin Bloom 
Assistant Professor of Aeronautical Engineering, Polytechnic 
Institute of Brooklyn, Brooklyn, N.Y. 


December 14, 1953 


— WRITER HAS PREVIOUSLY" ? COMMENTED upon certain as 
pects of the calculation of laminar boundary-layer stability 
with respect to small two-dimensional disturbances according to 
the Lin-Lees theory, citing an improvement in the procedure sug 
gested by Lees.’ Results similar to those shown in reference | 
have been independently obtained by Van Driest‘* with an anal 
ogous but slightly different procedure. In the course of pre 
paring for publication the work upon which references | and 2 
were based, it was found to be advisable to revise Lees’s pro- 
cedure still further. A brief statement of the revisions and some 
results, more completely given in reference 5, follow. (The nota- 
tion of Lees* is used here. Maximum, minimum, and critical 
values are denoted by the subscripts max., min., and crit., re 
spectively. ) 
Eqs. (19)* and (20)' are rearranged as follows: 
1 = ;'[1 — (26, — 1) + (4,2 + 0,;?)d\? X 
(1 +A)~!] 


u = &, (1 — (2 + 62 — 6,)¢,~'d] [1 — 2 (4 — 1d 4 
($;? + 2 — 2%, + 1)d7]~' (2) 


Gl Z, rd) (1 


where v = 2(c, Mo, T:), provided that a is not much greater than 
1 [cf. Eq. (24)]; alsoX = A(c, Mo, Ti). Plots of G versus z for 
various values of \ show that G > 0 may be double-valued in z 
when —0.1 < A < 7, that G < 0 may be double-valued when A > 
-(0.1, and that |G! has a minimum for these ranges of A (see 
Fig. 1). With v and ) calculated for given values of c, Mo, and 
7), the value or values of z which satisfy Eq. (1) can be found from 


the aforementioned plots of G. If no such value of z can be found 
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VALUES OF G REQUIRED TO CORRESPOND TO GIVEN VALUES OF . IF 


2 IS TO HAVE A SINGLE VALUE ; MINIMA OF IGI 


in a given case, no neutral disturbance that satisfies the necessary 
physical conditions exists in that case. Corresponding values of 
u follow from Eq. (2) and values of @ follow in principle from Eq. 
(27);3 the reliable calculation of a can be difficult if @ is much 
larger than 1. Appropriate values of the Reynolds Number R 
are derived from Eq. (26).4 

Of particular interest is the fact that, for given values of .1/o 
and 7}, the Cmaz. for a neutral subsonic disturbance generally 
occurs together with |G), in. rather than with %;, »,ar. (= 0.58) as 
assumed by Lees. However, for very small values of \ and for 
v>O, IG min. ANd ®;, »ar. Occur together, and Lees’s assumption 
is valid. This is important because of Lees’s further assumption 
that the minimum critical R (i.e., Rerie.) occurs together with 
Cmax. It is shown in reference 5 that C,,az, does not occur pre 
cisely with Rez. and that ¢,,,,, occurs simultaneously with a@»az 
only when A is near zero. (Lees’s curves, Fig. 4,° show that the 
values of R corresponding to a@»ar. are close to Reriz., but only 
small values of \ were encountered in the calculations he cited. 
Therefore, no positive correlation between Cnar, and Rep, could 
be established by the writer in reference 5 for general values of 
\ which may be fairly large. On the other hand, the absence 
of an approximate correlation between Car. and Rerir 


However, the solid curves in Fig. 1 


could not 
be demonstrated either. 
can be used to calculate the value of Car, for given values of 
Mo and 7;. Also shown in Fig. 1 is the curve suggested by 
Lees’ for calculating Cmar, with small values of \, and the curve 
previously suggested by the writer! for all values of \ but based 
upon the erroneous assumption that Car, and ®;, max, occur si 
multaneously for all values of X. 

It is also shown in reference 5 that the correlation between 
Rerit, ANd Cac. is not a factor in the calculation of surface tem- 
peratures that result in infinitely large values of Rerit.. This is 
because the necessary and sufficient condition for R-riz, to be in- 
1. The results of these 


finite is that ¢,,,,, Should approach 1 — .\/» 


calculations, which utilize the solid curves of Fig. 1, are shown 
in Fig. 2 

Only a tentative interpretation of the results in Fig. 2 can be 
offered for reasons explained in reference 5. In brief, it is con- 
sidered probable that flows in regions I, IV, and V are unstable 
at sufficiently high Reynolds Numbers, while flows in region II 
are probably stable. The flows in region III require further in 


t 


vestigation, although they have previously been considered t 


be stable.‘ Thus, the curves bounding regions IV and V 
interpreted as being possibly narrow bands of greatly increased 
stability. 

It has been authoritatively suggested to the writer that the 
curves bounding regions IV and V may indicate stability with 
respect to an additional mode of disturbance—that is, that an 
unstable mode might also exist for the indicated combinations 
of 7; and Mp. 
with the present form of the theory, since large values of ¢ a 
idditional u 


This possibility cannot be proved or disproved 








a would probably occur simultaneously for any 
stable mode in these narrow regions and therefore could not be 
treated with the present form of the theory. In this connection 
it can also be pointed out that Lees’s extension of the Heisenberg 
criterion® appears to be valid for large values of c, provided that 


a@ is very small, but may not be valid when both ¢ and @ are largt 


simultaneously 

In the work described above no attempt was made to consider 
the problem of improving the viscous solutions as discussed by 
Lin’ and by Cheng® or to consider the effect of the three-dimen- 


sional character of the disturbances 
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The Interaction of Two Similarly Facing Shock 
Waves 


G. V. Bull,* L. R. Fowell, and D. H. Henshawt 


Research Assistants, Institute of Aerophysics, University of Toronto, 
Toronto, Ontario, Canada 


December 10, 1953 


SUMMARY 


When one shock wave overtakes another traveling in the same direction 
in a medium for which y ;, the configuration resulting from the inter 
action consists of a reflected rarefaction wave, a transmitted shock wave 
and a contact surface separating the states behind these waves 

For y > 5/3 one may obtain the above result or a reflected shock wave 
with the contact surface and transmitted shock 

The analytical relationships enabling the determination of the final wave 
strengths as a function of the initial conditions are given, and graphs are 
presented showing the strengths of the transmitted shock wave and reflected 


rarefaction wave for the case of 1.40 


NOTATION 


(R, S, m, r,s) subscripts denoting the state of the physical system 
concerned 

S the undisturbed region ahead of the first shock wave 

R the region behind the second (overtaking) shock wave 

m = the region separating the first and second shock waves 
before interaction 

? = the region between the reflected wave and the contact 
surface after interaction 

s the region between the transmitted wave and the con 
tact surface after interaction 

5, 2 a shock or rarefaction wave, respectively, such that 

<—y the fluid enters the wave from the right 

SLR waves such that the fluid enters from the left 

+ 

. = contact surface 

p = pressure 

p = density 

7 specific heat ratio 

u a velocity measured relative to axes fixed in the physical 
system 

1/P pressure ratio across the system PR/ PS 

PRm = PR/ Pm 

Pms = Pm/ ps 

Ps Ps/ PS 

P. Pr/PR 

a (y + 1)/(7 1) 

B (y — 1)/2 

* Now Chief Aerodynamicist, Canadian Armament Research and De 


velopment Establishment, Quebec, Canada 
Tt Now Research Officer, National Aeronautical Establishment, Ottawa, 
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Fic. 1. The initial states in the physical system 
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Fic. 3. Physical states after interaction yielding reflected shock 


(1) RESULTANT WAVE CONFIGURATION 


NITIALLY, ONE CONSIDERS two shock waves traveling in the 
I same direction (Fig. 1). Relative to the right-hand shock 
the flow in (m) is subsonic, while relative to the left-hand skock 
it is supersonic; hence, the shocks must collide. 

A qualitative analysis in the (p, «) plane shows that this inter 
action may yield a reflected rarefaction wave (Fig. 2), a reflected 
shock wave (Fig. 3), or a reflected Mach line; while a ‘“‘trans- 
mitted” shock wave continues in the direction of the original 
waves in every case. 

In the quantitative discussion a one-dimensional flow was con 
being constant throughout any one 
density relationships 


sidered in an ideal gas (4 
system considered). The 
across rarefaction and shock waves yielded the fact that a con 
tact surface lies between the reflected and transmitted waves in 


pressure and 


every nontrivial case. 

A criterion was developed which, by means of the initial shock 
strengths and y for a given system, determines whether the r 
flected wave is a rarefaction or shock wave for that system 

The velocity increment across the initial system was written 

Ps a = ti 1 
uére—- us = 
\ ps 1+ aPns 

j Pata * 
‘ ns ~— 1) 4+ (1 — PPus) ( 

' Vp a + PP xs) 
while that across the transmitted shock in the case of a reflected 
Mach line is 


(1.1 


Ie — 1) ‘ 
uy — Us = Ps ta (1 P) (1.2 
Nos VPUP + a) 
Letting A designate the ratio (up — us)/(u, — us) one sees that 


A > 1 implies a reflected shock, A = 1 implies a reflected Mach 


line, while A < 1 corresponds to a reflected rarefaction wave 





P./Ps 


Transmitted Shock Strength 


KP, 


Strength of Reflected Rorefaction 


Mm 


Se 


I rare 


hock 


ted 


1.2 


hat 





READERS’ FORUM 

































































































































10000 -~ 
© : i Fam = Pos 
* 8000 Pam > fms ins > Frm 
| 
“ . P= 0.0001 
2 6000 
: l 
: : Sal | 
x 4000 an 
° LS SS a SS a SS eS > > > Ce ay Sa SS Fe 6 et ae eae eee 
500 P= 0.002 + Pon = Pin 
vo 400 (E quotion 2.1) Fam > Pms < Pms > Pam 
CY) 
= Parameters: ( 
: 300 One 
i) = | y s |.4 Pa e & 
5 ee or ‘npn enegesgeen ommpreree me. ee ae 
5 = RM ‘ms eee 
10.0 P=O.1 + Ra 
lo.) 
95 ike 
9.0 l og Badu l =a ser l =e — eee 
0.0001 000! 0.0! 0.1 1.0 
Reciprocal of Incident Shock Strength fe. 
Fic. 4. Transmitted shock strength in terms of incident shock strength. 
a” Tr T T T rE T a ae. y ~ =. 
: Rm>'ms ms > Rm 
2 
2 (Equation 2.2) 
® 
. — Parameters: 
50 Rm= eo. 
a 
} 
$) 
& 
~ 
t) 
x 
pe 
° 
oi5 
a 
c 
: 
Ps 
1.0 1 L Pi —_—— os l 
0.000) 0.00! 0.0! 1.0 


Reciprocal of Incident Shock Strength 


FG 


5 


0.1 
Pam 


Strength of reflected rarefaction wave in terms of incident shock strength. 











212 JOURNAL OF THE 
The condition A = 1 was written 


Pas — 1) RP = 1)? (RP, s? l 2a] x 


| 4 } 
a’ — 2a* +a l — — _ rh oe 0 ¢23) 
RP ns RP ms? RP ns) 
where k = ?p,,/P»,s. Considered as a function of y, this rela- 
tion yields a critical y = yc as a solution of the cubic in a. For 


y > yc a reflected shock occurs, while y < yc implies a reflected 


1 


rarefaction wave. Consideration concerning the smallest pos- 


sible physically admissible critical y showed that a reflected 
shock wave could not be obtained for 7 < 5/; under the assump 


tions of this analysis 
(Il) WAvE STRENGTHS 


The transmitted shock strength (7?,) was obtained in terms of 


the initial shock strengths in the following form: 
ye aP (’ ' Pa)(! | Pas 
Wit + a?) lo + Pm Fhe + Pos 
-§ [ft — Pm b+ ala 
od, ee 
‘ a a (: taPns) (| + Pee) , 
\ 1 a+ Pus Fa + Pn J) 
l 


(PP,)F — 


0 (23 


1.4 with P as 
The strength of the reflected rarefaction wave (1/7, 


Fig. 4 shows a plot of this relationship for 4 
parameter. 
was found to satisfy 


) ‘ B(1 + aP rm ) ] 
ae. 
Wie + Pe.) PP aun + 1) (Pan Y 
| ne l( PP x,, 4 a)PpR 
y rie 


PPPs, a4 
1) eee? _ 1 ~ PP, 


Va + aP ry 


and is shown in Fig 5 


The relationship giving the strength of the reflected shock for 
shock 


cases with 5/, was also found, and the reflected 


strength was plotted for a hypothetical case with y = 2 


(111) AppLlicaATIONS 


In wave interaction experiments the Eqs. (2.1) and (2.2) and 
the graphs of Figs. 4 and 5 enable one to specify the strengths of 
the waves resulting from the overtaking of a given shock wave 
by another of known strength. 


It will be noted that, for a given total pressure ratio (1//), the 
strength of the transmitted wave is greater the more completely 
this pressure rise can be made to occur over only one of the initial 
waves. For initial waves of comparable strength the trans- 
mitted wave strength is roughly one-half of the product of the 
initial pressure ratios, if the total pressure ratio is in the range 
500 to 10,000, and approaches this product in magnitude as 


(1/P) approaches one. 


An attempt to produce very strong shocks in a shock tube by 
this method, using two fairly strong initial shocks, has recently 
been brought to the authors’ attention. In this connection it is 
noted that for comparable initial waves a stronger transmitted 
wave should result when the first wave exceeds the overtaking 
wave in strength, because of asymmetry in this sense with re 


= Pns 


spect tothe line ? p, 
Experimental results have recently been obtained at the Insti- 
tute of Aerophysics confirming this analysis and will be discussed 


in a forthcoming report. 
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Remarks on Nonlinear Lift and Vortex 
Separation 


Hsien K. Cheng 
Aerodynamic Section, Bell Aircraft Corporation, Buffalo, N.Y 
December 1, 1953 


IT A RECENT NOTE, Adams! has given an interpretation of th, 
results of a theoretical investigation on vortex separation 
from the leading edges of slender delta wings,’ in which a simpli 
fied physical picture has been considered in that the free vor. 
tices separated from both edges are replaced by a pair of con 
centrated vortices in the cross-flow plane. Adams’ conely 
sion shows that, with the occurrence of the ‘free vortex pair,’ 
the result of the Jones low aspect ratio theory remains valid to 
the first order of approximation and that the nonlinear lift du 
to the vortex separation appears of the order of (A.R.)/ (a 
Nevertheless, his solution does not lead to a consistent deter 
mination of the nonlinear lift. Meanwhile, Edwards? points out 
that the inconsistency in Adams'’s results arises actually from ai 
error in the compatibility condition of the free vortices with their 
surrounding stream; by applying the correct condition, the sys 
tem of the vortex pair (together with their ‘‘feedlines”’) becomes 
dynamically free, and the inconsistency as to the question of in 
cluding the forces on the vortex feedlines is removed 

By using essentially the same model in describing the vortex 
separation phenomenon near the edge, the writer has obtained 
the litt of a rectangular plate at supersonic speed, which includes 
the contribution of the vortex separation from the side edge; 
in a similar manner, one can also determine, to the first order of 
a flat plate, taking 
the latter 


approximation, the impulse of motion of 
into account the phenomena of vortex separation ;‘ 
problem can again be interpreted through the cross-flow concept 
as one involving the “‘side-edge separation”’ from a low aspect ratio 
airfoil.6 Since results of investigation of this type are not gen 
erally known for the side-edge problem, it may be of interest to 
give a close examination of their characteristic features and t 
compare them with those that occurred in the other case pri 


viously mentioned. 


(1) THE VorRTEX PatR BEHIND A MOVING PLATE 


It will serve the present purpose to consider only an exampl 
of reference 4, in that a plate of width 2a submerged in an in 
finite expanse of fluid is set into transverse motion instantane 
1. The 


classical solution to the problem which admits infinite flow ve 


ously with a uniform velocity aU’, as illustrated in Fig 


locity and stream curvature at each of the edges is well known 
By introducing into the flow field a vortex pair, the infinite ve 
locities at the edges can be eliminated, and the solution is dé 
termined finally by applying the compatibility condition of the 
free vortex pair with its surrounding stream. Dynamical con 
sideration of the free vortices leads to the condition that the mean 
velocity of the group of vortices which has been shed before the 
instant ¢ be equal to the mean velocity of the surrounding stream 
at the same instant. An essential point contributed by Ed 
wards is that, in calculating the velocity of the free vortices, the 
path of the vortex center of the free vortex group under con 
sideration should not be confused with the locus of the c.g. of the 
total vortices in all time.* Let v denote the site of one of the 
“vortex pair” and T the total circulation around the free vortex 
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region; the motion of the vortex center of the free vortex group 
at the instant / may then be determined by 


dv* lj d dv 
= (lv) —a (1) 
dt I} dt dt 


The compatibility condition now reads as 


where the function on the right-hand side is the conjugate of the 
Both 
(2) and 


complex ‘“‘vortex-velocity” at the vortex center ¢ = v 
the circulation and the vortex-velocity appeared in Eq. 
are functions of y anda. The differential equation, Eq. (2), to 
gether with the initial condition v = a at ¢ = 0, is sufficient in 
principle to determine the location of the vortex pair as function 
of time and leads finally to the solution desired 

Accordingly, the equivalent solution for the delta wing prob 
lem of references 1, 2, and 3 appears to differ from the general 
case only at the last step—namely Eq. (2)—in that the first 
order differential equation becomes immediately integrable be 
cause of conical symmetry in ¢ (in particular, the width 2a being 
in that case a linear function of time, resulting in an expanding 
plate in a uniform stream). 

For the present case, the width 2a remains constant through 
out, and the integral solution of Eq. (2) is again seen possible 
by a separation of the variables. Moreover, if only the solution 
corresponding to small values of |v — a! is of primary interest, 
which implies a sufficiently small value of a, an explicit solution 
of the differential equation can be obtained by solving for « = 
¥ — a, which is the location of the free vortex center with respect 
tothe edge. Asa result of the boundary conditions given and 
the assumptions mentioned previously, Eq. (2) can be written, 


toits first order, as 
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ld 1 ia l Ve 
(le) = — tal . ( - 3 
I dt ‘ NB s/c" aa a? cad Ve — 


where 


r= al Va/2 ly eVe (v e+yv «| 

which yields the first order solution for € as 

e/a +i(talt/4a)y (4) 
where the negative signs correspond to the upward motion of 
the plate. It is seen from the above that the vortex pair trails 
directly behind the edges according to a time law of t’*. The 
corresponding value of the circulation around each of the free 
vortex pair, which measures its total strength, is 

I = 24Ua(alt/4a)'/*a (5) 
The solution is thereby completely determined. By virtue of 


the symmetry in the flow pattern, the total momentum received 
by the fluid medium may be evaluated simply as 


af x — 2at\ */+ 
rpa?a U + rpala? (6) 


a 


I(t) = 2xp lim i¢?W’(¢) 


C\— a 


AN INTERPRETATION OF THE ABOVE SOLUTION 


Obviously, the above solution pictures the cross-flow field of 
a slender airfoil in the event of vortex separation from its side 
edges. By multiplying the right-hand side by U and putting 
x/U for t, Eq. (6) provides a formula for the lift of a slender air 
foil possessing sharp side edges, which is correct up to the order 


° ro ro 4 Zax » J 
rpa?al? + relate (7) 
a 


where x is the Cartesian coordinate running downstream along 
the chord and its origin is chosen at the section where the vor 


considered—namely, 


lift |} 


tex separation first occurs 

Unfortunately, the present result interpreted from the previous 
solution implies the assumption of slenderness of the airfoil 
throughout the whole length of x, which will certainly limit its 
application to a great extent. However, in approaching the ex- 
treme of “infinitely small aspect ratio,’’ Eq. (7) must hold valid 
within the range of sufficiently small angle of attack a. 
upon comparison with a set of early data obtained by Winter,* 


the agreement appears to be almost perfect, in that an extremely 


Indeed, 


narrow rectangular flat plate of aspect ratio '/3) was used (see 
Fig. 1). For rectangular plan form of higher aspect ratio, say 
between 1/2 and 1/10, the discrepancy becomes appreciable 
(the data for aspect ratio of 1/4 appeared in the figure are taken 
from reference 5), yet Eq. (7) still accounts for the total nonlinear 
effect to a substantial degree; 
has also been qualitatively confirmed by Winter’s measurement 


the corresponding lift distribution 


Furthermore, within the degree of accuracy afforded by the pres 
ent theory, the equation of the free vortex core in space deduced 
directly from Eq. (4) seems to agree also with a side-view photo 
graph in reference 6 

The discrepancy between theory and tests that occur in the 
is perhaps due mainly to the separa 
and the 


higher aspect ratio range 
tion phenomena at both the vicinities of the leading edge 
trailing edge of the rectangular plate, which are in fact the real 
stalling phenomena that cause a decrease in lift. Since for a 
narrow rectangular plate the nonlinear effect contributed by the 
latter phenomenon is not large in comparison with the effect due 
to side-edge separations, an empirical formula may be adopted 


after Eq. (7) for the higher aspect ratio range as 
C, = (#/2) (A.R.)a + ka? (8) 
where the coefficient of the linear terms may be modified 


cording to reference 5 for aspect ratio higher than 1/2. It is 
interesting to note that the nonlinear lift coefficient deduced 
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from Eq. (7), as presented in the figure, happens to differ from 
either formulas of Adams and Edwards by a factor of (2)'/? 


(3) COMPARISON WITH THE DELTA WING SOLUTION 


In view of the fact that the ordinary solution of the flow around 
a leading edge possesses a singularity of a higher order of in- 
finity than around a side edge, it is to be anticipated that a 
stronger separation phenomenon may occur in the event of the 
leading-edge separation. Indeed, this is revealed by comparing 
the solutions of the separation models themselves, for, according 
to the solution above, the total strength of the free vortex is of 
the order a‘/* with a typical dimension of the region in the order 
of a’/*, whereas close examination of the solution of references 
1, 2, and 3 shows that the strength is of the order a and the size 
of the vortex region of order a’/*. To be sure, for the narrow 
delta wing problem, 


P = rUaa (9) 


€ 1 a ae (2 
-_— 2 \4da ; da 


Inasmuch as the nonlinear lifts belong to the same order of 
magnitude in both cases—that is, the order of (A.R.)'/%a’ 
whereas the strengths of the vortex region differ, it seems possible 
that the error implied in Eq. (7) could be less than that which 
occurred in the corresponding formula for the delta wing. 

Finally, it should be noted that the power law ¢/* or x’/? for 
the vortex core has likewise occurred in the case of the rolling up 
of the trailing vortex sheet behind a lifting surface.’ Further 
remarks and interpretations of the results of the similar nature 


(10) 


3 


can be found in reference 4. 
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Remarks on ‘‘Subsonic Compressibility 
Corrections for Propellers and Helicopter 
Rotors’’! 


M. Z. v. Krzywoblocki 
Professor of Gasdynamics and Theoretical Aerodynamics, 
University of Illinois, Urbana, Ill. 


December 3, 1953 


iy REFERENCE | THE AUTHORS came to the conclusion that the 
assumptions of previous authors” * concerning the axial inter- 


ference velocity factor due to induced axial flow are incorrect 
They object to the assumption made in reference 2 that p: = p, 
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(no change in density through the propeller disc), on the basis 
that for a perfect gas any pressure discontinuity in an adiabatic 
flow must be accompanied by a corresponding density discon 
tinuity. They further object to the assumption made in refer 
ence 3 that p; = po, where p; refers to slipstream conditions at 
an infinite distance downstream, on the basis that one cannot as 
sume the existence of a thin actuator disc. As a possible soly 
tion, they propose to use the exact incompressible flow relations 
as the first step of an iterative process; but, they neither go be 
yond that first approximation nor discuss the iteration process 
nor mention anything about it. Hence, the following facts mus 
be stated: 

(i) The first-order approximation proposed by the author re 
fers fundamentally to the subsonic incompressible flow and to 
that region of compressible flow where one can neglect higher 


order terms. 


(ii) The above implies that the authors did not find a solution 
of problem in question in the region where the compressibility 
phenomena play a full part and where the first-order term is not 
strong enough; they object to previous suggestions but do not 
offer any improvement. 

(iii) Concerning the iteration process, there may be an equal 
probability of obtaining equivalent results by assuming that p, 
is a first approximation to p, (reference 2) or that po is a first ap 
proximation to p; (reference 3) and performing the iteration proc- 
ess for p2 or p;. There is absolutely no obvious reason why the 
iteration process for 6; in terms of 6; proposed by the authors 
should be superior to the iteration processes for p» or p3. 


Now a few words about the iteration process. Any iteration 
process may give acceptable results, but very often the domain 
of the validity of that process is small. Thus, it seems possible 
that the iteration process on p2 in terms of p;, limited to the thin 
actuator disc, may give acceptable results. Similarly, the process 
for p; in terms of po, limited to some narrow domain, may be ac 
In the past the author of the present remarks applied 


ie., they di 


ceptable. 
a few iteration processes with a negative result 
verged. It looks as if, in the problem in question, there were 
some inherent properties repudiating even the law of an iterativ¢ 
process and may be subject only to some experimental pro 
cedure. This is a situation similar to a system having more 
unknowns than available equations, and here may be an ex 
planation of the failure of iterative processes. But, of course, it 
does not mean that some way will not be found to solve the 


problem analytically. 
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On the Curved Shock Wave Due to an Infinite 
Wedge 


Kinzo Hida 
Naniwa University, Sakai Near Osaka, Japan 
December 14, 1953 


Ss’ FAR IT HAS BEEN KNOWN that there are two possible com 
figurations of the shock wave due to an infinite wedge placed 
in a supersonic uniform flow, provided that its vertex angle (260) !s 


smaller than a certain critical value. Both are straight shock 
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waves and are called ‘‘strong-shock”’ and ‘‘weak-shock,"’ respec- 
tively. 

In this paper we shall show, with the aid of analysis in the 
hodograph plane, that two more shock configurations are pos 
sible. They are both curved shocks; one is an attached shock 
that is weak at the vertex and becomes strong at infinity, while 
the other is a detached shock changing from normal-shock to 
strong one at infinity. 

When the Mach Number \/, of the uniform flow at infinity 
upstream is prescribed, not only the magnitude and the direction 
of the velocity but also the initial curvatures of stream lines im- 
mediately behind the shock are given by a shock polar as well as 
the shock spines (see Fig. 1).? 

Let a straight line be drawn through the origin 0, making an 
angle # with OP and let A, B denote the points of intersection of 
such a straight line with the shock polar. Then these two points 
correspond to the weak-shock and strong-shock, respectively 
In these cases the flow behind the shock wave becomes uniform 
again, and the flow field shrinks into a point in the hodograph 
plane. But we can show two other possibilities of shock con 
figuration by considering the solutions that correspond io the re- 
gion BON (II) or the region BMA (1). 

For brevity’s sake, we take the Laplace equation Ay = 0 as 
the fundamental equation for the stream function ¥ in each re 
gion. (This simplification is justified if use is made of a suitable 
afine transformation or the W.K.B. method suggested by Prof 
|, Imai.?) 

Boundary conditions are as follows: 


\on AB (for the region I), 
jon BO and ON (for the region I1). 

(ii) Initial direction of every stream line should coincide with 
that of the spine at each point on the shock polar. 


i)v = 0 


To find appropriate solutions we shall discuss the two cases 


separately 


(i) REGION I (WEAK-STRONG SHOCK) 


In this region we shall find the solution as shown in Fig. 2. 


Considering a corner flow and a generalized doublet flow in the 
close vicinity of A and B, respectively, and superposing them, we 


get 





Fic. 1. 











Fic, 2. (Ordinates are 


Reduced hodograph plane for region I. 
stretched twice as large as abscissae. ) 
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Reduced hodograph plane for region II. (Ordinates are 
stretched twice as large as abscissae. ) 


Fic. 3 





(4) (2)(3) 





Fic. 4. Possible shock configurations due to an infinite wedge. 
(1) Weak shock; (2) strong shock; (3) weak-strong shock; (4) 
normal-strong shock. 
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Fic. 5. Probable flow patterns tending to our cases as limit 


VW = Cr sin m6 + Cop" sin n(x 3) 


where (7, @) and (p, #) are the polar coordinates with the initial 
line BA and the origin at A and B, respectively. If indices 
m, n be determined so as to satisfy the second boundary condition 


(ii) at A and B, we have 
nu = Br bp - &R nm = | Ba 54 


Moreover, we can determine the ratio of arbitrary constants, C,+ 
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At first sight these flow patterns seem to be rather curious 
because of having a certain characteristic length (e¢.g., the radiys 
of curvature of the shock wave at its nose), notwithstanding that 
a ee a ee a = infinite wedge has no such ve _ may, however, 

ye interpreted at least qualitatively as follows: 
mpy'™ sin (Bu =e og 0x} Taking the solution with the sink S at a certain point on BA or 
BO as shown in Fig. 5 (i), (ii), the corresponding flow patterns 
(ii) REGION II (NORMAL-STRONG SHOCK ) may be expected to be those given in Fig. 5 (iii), (iv). These flow 
patterns will be reduced to either of the above-mentioned cases 


C2, in such a way that the boundary condition (ii) is satisfied at 
one more point, say ./ (i.e., point of maximum deflection). We 


then have 


To obtain the solution having the singularity at B as shown in 
by making S approach B. 


The detailed discussions with some numerical calculations 

v = Im[C/(s™ + b™)"| will be reported elsewhere in the near future. 
In conclusion, the author wishes to express his cordial thanks to 
Prof. S. Tomotika, of the University of Kyoto, and to Dr, K. 


m = 1/80, n = (Bp 6p)/(m + bn 9) Tamada for their continual guidance and encouragement through- 


Fig. 3, we take 


with 


r , . ) his work. 
Here, 2 = ge’? is the complex variable in the hodograph plane and out this work 
b denotes the magnitude of 0B. This function takes the value REFERENCES 


Wv = 0on NO and OB and satisfies the condition of the spine at B. 
' Busemann, A., A Review of Analytical Methods for the Treatment of Flows 


(The second condition (ii) at N is automatically satisfied. ) 

aie ( : ‘ with Detached Shocks, NACA TN No. 1858, 1949. 

Fig. 4 shows an example of four possible shock configurations 2 Imai, I., Application of the W.K.B. Method to the Flow of a Compressitle 
0 and 0) = 11° 45’. Fluid, I, Jour. Math. Phys., Vol. 28, pp. 173-182, 1949. 


for the case when MM, = 
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